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1.  Objectives 

The  basic  goal  of  our  research  effort  has  been  the  development  of  computational  methods 
and  tools  which  optimally  exploit  the  analytical  procedures  natural  to  aerodynamic  theory.  This  has 
resulted  in  a  variety  of  procedures  of  non-standard  form  for  treating  a  wide  range  of  problems  in 
gas  dynamics.  We  believe  that  our  research  effort  has  made  significant  advances  in  subsonic, 
transonic  and  supersonic  gas  dynamics. 

Specific  objectives  of  the  program  have  been;  the  use  of  natural  coordinates,  e.g.  streamlines, 
characteristics,  "potential"  lines  and  so  forth  in  the  formatting  of  compressible  computer  codes; 
solution  of  general  inverse  or  design  problems  in  aerodynamics;  use  of  machine  algebra  to  format 
codes  and  deal  with  non-standard  problems;  the  development  of  a  method  of  parametric 
differentiation  to  extend  generally  existing  codes  (in  addition  to  ours)  to  continuous  ranges  of 
validity  in  parameter  space  (i.e.  Mach  number,  thickness  ratio,  camber  and  more  general 
parameters  specifying  a  body  shape.) 

All  of  the  stated  objectives  in  the  original  three  year  proposal  have  been  accomplished.  In 
addition,  as  described  in  the  following  section,  we  have  also  achieved  a  number  of  extensions  and 
additional  results  not  anticipated  in  the  original  proposal. 


2.  Research  Narrative 


Our  research  effort  began  with  a  method  for  treating  two-dimensional  supersonic  flows  past 
airfoils.  This  is  based  on  transformation  to  streamline  and  principal  characteristic  coordinates,  and 

results  in  a  rapidly  convergent  and  accurate  solution.  Both  body-fit  and  shock-fit  coordinates  are 

generated  by  this  method.  The  codes  produced  by  our  method  are  perhaps  the  most  efficient 
now  in  existence.  A  typical  calculation  takes  a  small  fraction  of  a  second  on  a  mainframe, 

[reference  1] 

This  method  was  next  extended  to  the  inverse  or  design  problem  for  two-dimensional 

supersonic  flows.  Through  analytical  procedures  the  inverse  problem  was  transformed  to  a  direct 
problem  of  different  type.  The  result  is  a  speedy  accurate  procedure  for  determining  shape  from 
a  given  pressure  distribution,  [reference  2] 

A  method  based  on  streamlines,  characteristics  and  Riemann  functions  has  also  been 

introduced  for  supersonic  flow  over  axisymmetric  bodies.  Starting  from  a  simple  approximation,  an 
iterative  procedure  is  developed  which  converges  rapidly  to  the  exact  solution.  The  scheme  is 

both  body-fit  and  shock-fit.  As  a  result,  the  procedure  is  computationally  efficient,  inherently 

accurate,  and  requires  relatively  few  points  to  calculate  the  entire  flow  field.  Both  the  direct  and 
inverse  design  problems  are  treated.  For  a  thin  axisymmetric  body  traveling  at  low  supersonic 
Mach  number,  our  results  show  the  presence  of  a  pressure  minimum  on  the  body,  a  phenomenon 
which  seems  to  have  gone  unnoticed,  [reference  3] 

The  method  described  for  the  axisymmetric  case,  has  been  adapted  to  the  treatment  of 
non-axisymmetric  bodies.  In  particular,  we  consider  flow  in  azimuthal  planes  and  develop  a 
procedure  based  on  near  characteristics  and  projected  streamlines.  The  cross-talk  between 

azimuthal  phases  defines  the  basis  of  an  iteration  procedure  which  is  rapidly  convergent.  As  is 
the  case  for  axisymmetric  flow  relatively  few  computational  points  are  required.  Both  the  direct 
and  indirect  flow  problems  have  been  treated,  [reference  4] 

We  have  treated  subsonic  gas  dynamics  in  the  tangent  gas  approximation.  Using  a  highly 
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analytic  basis  a  very  fast  and  accurate  method  of  solution  has  been  developed  for  the  numerical 
solution  of  subsonic  problem.  Comparison  of  tangent  gas  and  exact  flows  show  that  the  former  is 
extremely  accurate  except  at  locations  that  are  critical.  Tangent  gas  solutions  when  used  as  the 

first  step  in  the  iterative  solution  of  the  excat  flow  field  are  shown  to  give  substantial  reduction 

in  computation  time,  [reference  5] 

The  inverse  problem  in  the  tangent  gas  approximation  has  also  been  considered,  and  an 
exact  method  for  designing  airfoils  developed.  Constraints  on  the  speed  distribution  are  easily 

implemented.  A  simple  numerical  algorithm  which  is  fast  and  accurate  has  been  obtained. 

Comparison  of  designed  airfoils  using  the  tangent  gas  method  with  exact  Euler  results  is  found  to 
be  excellent  for  subcritical  flows,  [reference  6] 

The  methods  used  in  the  treatment  of  the  tangent  gas  have  been  extended  to  the  full 
two-dimensional  potential  equations.  A  powerful  combined  analytical  and  numerical  procedure  now 
permits  the  treatment  of  both  the  direct  and  inverse  problem  for  subsonic  and  transonic  problems. 

This  method,  which  still  needs  further  implimentation,  may  have  a  significant  impact  on  the  way 

that  transonic  airfoils  are  designed  in  the  future,  [references  7  and  8] 

The  flow  of  an  inviscid,  irrotational  and  compressible  perfect  gas  in  the  upper  half  plane  is 
used  as  a  model  to  investigate  the  transonic  controversy.  The  solution  of  the  complete  potential 
equation  for  the  velocity  potential  tf>(x,y),  with  boundary  condition:  <p  +  c  <Py  =  U  sin  x  on  y  =  0, 

is  developed  as  a  regular  perturbation  series.  36  terms  of  the  series  are  determined  by  computer. 

The  effective  boundary  condition  is  varied  with  the  choice  of  c;  and  for  each  of  the  velocity 
series,  its  nature  and  the  location  of  the  singularity  nearest  to  the  origin  are  investigated  using 

the  ratio  method  of  Domb  and  Sykes  and  Pade  approximants.  The  result  of  the  analysis  shows 

that  the  phenomenon  of  shockless  transonic  flow  is  dependent  on  the  imposed  boundary 

condition-which  for  this  example  is  the  constant  c.  The  relationship  of  series  convergence  to 
local  sonic  conditons  shows  no  obvious  pattern.  Cases  for  which  convergence  lies  below,  above  or 
is  at  critically  were  found.  Moreover,  the  connection  of  divergence  to  the  appearance  of  shocks 
is  also  not  apparent.  For  one  class  of  flows  divergent  series  could  be  resummed  to  yield 


shockless  conditions  for  all  Mach  numbers.  Significant  use  of  the  machine  algebra  code,  Macsyma, 
was  used  in  this  study,  [reference  9] 

We  have  also  treated  steady,  inviscid  supersonic  flow  over  three  dimensional  wing-like  bodies 
numerically  as  a  coupled  set  of  two-dimensional  characteristic  problems.  Shock  fitting  is  used  in  a 
boundary  fit  coordinate  system  and  the  calculation  is  second  order  accurate.  The  difference 
equations  are  solved  iteratively  and  the  use  of  an  accurate  approximation  step  results  in  rapid 
convergence.  A  variety  of  different  iteration  methods  are  considered  and  compared.  Incorporation 
of  a  flexible  data  structure  in  the  program  allows  for  an  efficient  use  of  memory  and  allows  a 

wide  range  of  wing  geometries  to  be  handled.  Results  for  tapered,  delta  and  swept  wings  at 
several  Mach  numbers  are  compared  with  two-dimensional  theory.  The  technique  is  applied  to 

both  the  direct  and  inverse  problems.  Derivations  are  carried  out  in  a  general  manner  allowing 
extensions  of  the  method,  [reference  10] 

We  have  developed  a  method  using  parametric  differentiation  which  can  significantly  extend 
any  numerical  study.  In  brief  flow  past  a  body  is  in  general  specified  by  a  variety  of  parameters 
such  as  thickness,  angle  of  attack,  camber,  Mach  number  as  well  as  others.  A  particular  flow  is, 
therefore,  characterized  by  a  single  point  in  the  corresponding  parameter  space.  Conversely,  the 

numerical  calculation  of  a  particular  flow  field  yields  information  at  just  one  point  of  the 

parameter  space.  However,  the  nature  of  a  continuous  range  of  nearby  flow  fields  is  of 

fundmental  significance  in  the  design  and  performance  of  aircraft.  To  treat  this  generally,  one  can 
consider  the  variational  equations  (which  are  linear)  obtained  by  differentiating  the  exact  equations 
with  respect  to  each  of  the  relevant  parameters.  The  resulting  matrix  of  derivatives  of  flow 
quantities  is  referred  to  as  the  Jacobi  matrix.  The  subsequent  procedure  is  in  principle  now 

straightforward.  One  integrates  the  nonlinear  governing  equations  --  which  results  in  the 
determination  of  just  one  point  in  parameter  space  --  and  simultaneously  the  variational  equations 
governing  the  Jacobi  matrix.  The  last  is  then  used  to  describe  (he  neighborhood  of  the  already 

determined  point  of  the  parameter  space.  Since  the  variational  equations  are  linear  the  additional 

computational  time  required  for  their  integration  is  modest. 


Frequently,  when  calculating  the  flow  about  a  body,  one  is  interested  in  how  the  flow  will  change 
if  the  base  configuration  is  altered.  For  example,  one  may  want  to  know  what  will  happen  at  a 
slightly  different  angle  of  attack,  wing  loading,  camber  or  thickness.  To  answer  such  questions 
each  parameter  change  is  traditionally  considered  as  a  separated  case  and  flow  simulation  code  is 
repeatedly  run.  It  could  be  argued,  quite  effectively,  that  in  many  instances  this  is  not  an 
efficient  use  of  resources.  Why  undertake  an  entirely  new  calculation  of  the  flow  when  we  know 
the  results  at  a  nearby  state?  The  method  which  we  have  developed  allows  efficient  generation 
of  solutions  in  the  neighborhood  of  a  base  solution. 

Thus  far  we  have  applied  the  Jacobi  matrix  technique  to  five  problems.  The  direct  calculation  of 
inviscid  supersonic  flow  about;  two  dimensional  airfoils  of  varying  thickness,  angle  of  attack  and 
camber;  axisymmetric  bodies  of  varying  thickness  and  taper:  and  the  design  (inverse)  calculation 
of  inviscid  supersonic  flow  past;  airfoils  described  by  a  given  family  of  pressure  distributions; 
axisymmetric  bodies  described  by  a  given  family  of  pressure  distributions.  Also  to  subsonic 
potential  flow  about  two  dimensional  airfoils  by  modifying  FL036.  Results  of  these  calculations 
show  that  Jacobi  method  allows  for  the  efficient  and  accurate  generation  of  parametric  solutions 
in  the  neighborhood  of  a  known  solution,  [references  11  and  12] 
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An  approximate  solution  is  developed  for  two-dimensional,  steady,  Im viscid  super¬ 
sonic  How  over  an  airfoil.  This  approximation  produces  accurate  results  for  a  wide 
range  of  Mach  numbers  and  airfoil  thicknesses.  It.  is  list'd  as  the  starting  point  for  a 
rapidly  convergent  iterative  numerical  solution  of  the  exact  equations.  A  co-ordinate 
system  consisting  of  the  principal  characteristics  and  streamlines  is  employed. 
Examples  computed  for  a  symmetric  airfoil  reveal  several  interesting  features  in  the 
tail  shock  and  the  flow  behind  the  airfoil. 


1.  Introduction 

In  this  paper  we  consider  the  computation  of  inviscid  supersonic  (low  over  a  two- 
dimensional  airfoil.  While  the  final  step  in  our  investigation  is  numerical,  we  attempt 
to  incorporate  as  much  as  possible  our  analytical  and  physical  knowledge  of  such 
flows.  The  approach  is  well  suited  both  for  numerical  integration  and  for  the  inter¬ 
pretation  of  the  resulting  flow  phenomena.  A  preliminary  version  of  this  approach 
for  the  case  of  one-dimensional  unsteady  flow  has  already  been  reported  (Sirovich  & 
Chong  1980;  Chong  &  Sirovich  1980).  In  the  present  investigation  several  new  or 
little-known  effects  concerning  the  tail  shock  and  flow  behind  a  two-dimensional  air¬ 
foil  emerge.  These  are  discussed  in  §6. 

There  are  two  main  nonlinear  approximations  for  this  problem.  Small -amplitude 
theory  gives  solutions  valid  provided  the  airfoil  thickness  is  not  too  great  and  the 
Mach  number  is  not  too  high.  Under  these  conditions  the  leading  shock  wave  is  fairly 
weak  and  the  solution  is  approximately  given  by  a  simple  wave  involving  only  the 
characteristics  emanating  from  the  airfoil  (Friedrichs  1948;  Light.hill  1900).  Variations 
in  the  entropy  and  in  the  Riemann  invariant  which  is  carried  along  the  down  running 
characteristics  are  only  of  third  order  in  the  shock  strength,  so  t  he  resulting  approxima¬ 
tion  is  valid  to  second  order.  A  correction  in  the  tail  shock  region  is  necessary  to  obtain 
a  second-order  solution  there  (Caughey  1909). 

The  second  tyjie  of  approximation,  shock  expansion  t  heory,  originated  by  Epstein 
(1931 ),  employs  the  fact  that  even  for  flows  with  strong  shocks,  for  which  the  assump¬ 
tions  of  small  perturbation  theory  do  not  hold,  the  effect  of  the  down  running  character¬ 
istics  remains  small.  This  leads  to  an  analytic  solut  ion  at  the  airfoil,  which  has  been 
generalized  by  several  authors  (Eggers,  Syvertson  Kraus  1953;  Meyer  1957)  to 
provide  approximate  solutions  for  the  entire  flow  field.  In  another  approach,  Jones 
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(1903)  lias  derived  by  a  perturbation  met  bod  an  approximate  solution  between  simple 
wave  theory  and  generalized  shock  expansion  theory. 

In  §  4  we  derive  an  approximate  solution  which  is  closely  related  to  these,  but  which 
applies  its  assumptions  more  consistently  and  is  somewhat  more  accurate.  This 
approximation  includes  both  shock  expansion  theory  and  the  second -order  theories 
of  Friedrichs  and  Caiighev.  'flic  derivation  and  the  numerical  computation  of  the 
solution  are  facilitated  by  the  use  of  the  principal  characteristics  and  the  streamlines 
as  co-ordinates  (§  3).  Adamson  (1908)  lias  used  a  similar  co-ordinate  system  in  another 
context.  For  a  problem  in  which  the  down-running  characteristics  are  also  important 
(e.g.  flow  in  a  nozzle),  this  approach  is  less  appropriate. 

The  approximate  solution  is  used  as  the  starting  point  for  an  iterative  numerical 
computation  of  the  exact  solution  (§  5).  The  high  accuracy  of  the  approximation  leads 
to  the  exact  solution  after  only  a  few  iterations.  This  procedure  is  different  from  most 
numerical  methods  for  hyperbolic  problems.  Typical  methods  apply  one  of  a  variety 
of  differencing  schemes  (for  a  comparison  of  several  such  schemes  see  Taylor,  Ndefo  & 
Masson  1972)  to  the  equations  in  their  standard  form  and  compute  the  solution  by 
‘marching’  along  in  the  downstream  direction.  One  disadvantage  of  these  methods  is 
that  at  low  Mach  numbers  short  step  sizes  are  required  for  stability.  The  method  of 
characteristics  (Liepmann  &  Roshko  1957,  eha.  12)  can  also  be  used  for  this  problem, 
although  it  is  considered  in  general  to  be  somewhat  unwieldy  for  machine  computation. 
The  BVLR  method  (Babenko  et  al.  19Gfi:  Holt  1977)  is  a  finite-difference  method 
which  is  partly  based  upon  the  method  of  characteristics.  The  transformation  of  co¬ 
ordinates  employed  here  also  results  in  a  method  which  is  closely  related  to  t  he  method 
of  characteristics. 

Special  account  must  be  taken  of  the  appearance  of  shock  waves  in  this  type  of 
problem.  In  finite-difference  methods  this  can  be  done  through  shock-capturing 
difference  schemes,  or  through  explicit  shock  fitting  (e.g.  Salas  1970).  In  the  present 
method  the  shock  waves  can  be  naturally  incorporated  in  the  new  co-ordinate  system 
as  fixed  boundaries  of  the  flow  field. 

2.  Formulation  of  problem 

We  consider  uniform  flow  of  Mach  number  Jl/0  >  1  incident  upon  a  two-dimensional 
airfoil  (see  figure  1).  It  is  assumed  that  there  arc  attached  shocks  at  the  leading  and 
trailing  edges,  and  that  the  flow  remains  supersonic  everywhere.  The  flow  fields  above 
and  below  the  airfoil  can  be  computed  independently,  up  to  the  appearance  of  the 
tail  shocks.  The  tail  shock  and  the  flow  behind  it  for  the  ease  of  a  symmetric  airfoil 
are  treated  in  appendix  B. 

The  co-ordinates  x  and  y  arc  scaled  by  the  airfoil  length;  the  pressure  p  and  the 
density  p  by  their  upstream  values  pa  and  p0;  the.  velocity  (it,  v)  =  {q  cos  0,  q sin  0)  and 
the  speed  of  sound  a,  by  the  upstream  speed  of  sound  n0;  and  the  entropy  s,  which  is 
set  to  zero  upstream,  by  the  gas  constant  R.  We  consider  a  perfect  gas  with  constant 
specific  heats  c,.  =  li/{y  —  1)  and  cp  =  ye,.,  for  which  the  equation  of  state  is 

V  =  f,y cxp | (y  -  l).s] 

and  the.  speed  of  sound  is  given  by  a-  =  p/p.  'flic  calculations  here  were  done  for 
y  —  1-4.  Modifications  for  the  case  of  a  gas  with  a  general  equation  of  state  are  out¬ 
lined  in  appendix  A. 
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Figure  1.  Supersonic  flow  over  a  symmetric  25%  thick  circular  arc  airfoil 
at  upstream  Mach  number  Mu  =  2-5. 


The  equations  of  inviscid  two-dimensional  steady  flow  are  conveniently  written 
with  the  entropy  s,  the  flow  angle  0,  and  the  Mach  angle  fi  =  sin_1(l/A/)  (where 
M  =  q/a  is  the  local  Mach  number)  as  dependent  variables.  All  other  physical  quan¬ 
tities  can  be  obtained  from  these  and  Bernoulli’s  equation 


«2+^?2= i+v^- 


The  equations  of  motion  in  characteristic  form  are  (Meyer  I960,  p.  273) 


ds  =  0  on  streamlines  =  tan  0 ; 

dx 


d(0  ±  PUi))  =  ±  S'“-2—  do  on  C±  ^  =  tan  ( 0  ±  /t) ; 
2y  ax 


where  P(/i)  is  given  by 

P(p)  =  AHan_I(AHan/t) -p,  A  =  (y  +  1  )/(y  —  I). 


The  streamlines  and  the  C+  characteristics  are  shown  in  figure  1.  The  quantities 
r±  =  0  ±  P(/i)  are  called  the  Rieinann  invariants. 

If  the  airfoil  surface  is  specified  as  y  =  f{x),  the  appropriate  boundary  condition 
there  is 

tan#  =  f'(x)  on  y  =  /(.r).  (4) 
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The  jumps  in  0,  /i  and  s  across  a  shock  are  governed  by  the  Rankine-Hugoniot  con¬ 
ditions  (Liepmann  &  Roshko  1957,  p.  85).  All  three  quantities  can  be  written  as 
explicit  functions  of  il/0,  y  and  the  shock  angle,  tj. 

3.  New  co-ordinate  system 

As  mentioned  in  the  introduction,  in  a  problem  with  weak  shock  waves  deviations 
in  s  and  r~  from  their  upstream  values  are  third-order  quantities.  This  is  shown  in 
figure  2,  where  As  and  A r~  are  plotted  on  a  logarithmic  scale  against  the  deflection 
angle  0,  for  several  Mach  numbers.  As  0  -*■  0,  the  curves  approach  straight  lines 
of  slope  3.  While  As  and  Ar~  are  both  third-order  quantities,  for  a  given  Mach 
number  the  jump  in  r~  is  always  significantly  smaller  than  that  in  s.  This  suggests 
that  for  weak  to  moderate  strength  shock  waves  the  flow  field  can  be  considered 
primarily  an  interaction  between  a  simple  wave  and  an  entropy  variation,  with  r~ 
playing  only  a  small  role. 

This  leads  us  to  introduce  a  co-ordinate  system  (a,/?)  consisting  of  the  streamlines, 
a  =  constant,  and  the  principal  (C+)  characteristics,  /?  =  constant.  Taking  a  and  fi  as 
the  independent  variables,  x  and  y  must  satisfy 

Up  =  xp  tan  0,  yx  =  .ra  tan  (<?+//) 

The  entropy  equation  (2)  becomes 

sfi  =  0, 

or  s  =  s(a).  Equations  (3+  )  and  (3—  )  become 

n,  sin  2/i 
(0  +  P(/i))x  =  -—-s' {a) 

and 

(5+-S)  «'-'«> 

where 

Xp 

Using  (7),  equation  (8)  can  be  simplified  to 

(0- P(ji))p  =  (I  -  tan  Otan/i)  ~£  0  (10) 

Equations  (5)— (7)  and  (10)  are  five  equations  in  five  unknowns:  0,  ft,  s,  x  and  y. 

The  boundary  and  shock  conditions  in  the  a.p  plane  can  be  simplified  by  normalizing 
a  and  /?  appropriately.  We  let  the  airfoil  surface  be  the  st  reamline  a  =  0,  and  normalize 
P  by  setting  /?  =  x  at  a  =  0.  The  boundary  condition  (4)  then  becomes 

x(0J)  =  p,  y(0,P)  =f(P),  0(0,P)  =  tan-«/'(/?).  (11) 

One  convenient  way  of  normalizing  a  is  to  take  the  front  shock  a  ngle  ?;(a)  to  be  given  by 

tan ^(a)  =  (1 -a)  tan  »/(0)  +  a  tan//0,  (12) 

where  7/(0)  is  known  from  solving  the  shock  conditions  at  the  leading  edge,  and  /i0  is 
the  upstream  Mach  angle,  which  the  shock  approaches  far  away  from  the  airfoil. 


1  —  tan  0  tan  /t ' 


(5) 

(6) 
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Fiqurk  2.  Jumps  in  entropy  s  and  Riemann  invariant  r~  across  a  shock  wave  as  functions 
of  deflection  angle  0,  at  various  Mach  numbers: - ,  A*; - ,  A r~. 

Hence  a  =  1  corresponds  to  x,y->co.  If  y  is  not  a  strictly  decreasing  function,  a 
different  normalization  must  be  used.  The  flow  field  in  the  upper  half-plane  is  mapped 
into  a  finite  region  in  the  a/?  plane,  as  shown  in  figure  3.  The  principal  characteristics 
become  vertical  lines,  and  the  streamlines  become  horizontal  lines.  The  front  shock 
maps  into  some  curve  /9(a),  and  the  left-  and  right-hand  sides  of  the  tail  shock  into 
two  separate  curves  /?2(a)  and  /?3{a).  The  discussion  of  the  tail  shock  is  left  to  appendix 
B.  With  the  shock  angle  y(a)  agiven  function,  the  shock  conditions  can  be  immediately 
solved  for  0 (a, /1(a)),  /r(a, /?(<*)),  and  s(a).  The  shock  /9(a)  itself  will  in  general  depend 
on  the  rest  of  the  solution,  however. 

It  is  possible  to  eliminate  y  from  the  equations  by  setting  yafl  =  in  (5).  Using 
(10),  this  gives 

0  =  xafl/^  +  i/l  +  ^(/l))flc  ot/i  +  (0  +  /i)fit&n(d  +  /i),  (13) 

which  can  be  integrated  to 

.r(a,/?)  =  x(0,/?)  +  j*  A{a) a~x  cos  (0+/i) da,  (14) 
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Figure  3.  Flow  field  corresponding  to  figure  1  in  a/1  plane.  Streamlines  map  into  horizontal 
lines,  <x  =  const.,  and  C +  characteristics  into  vertical  lines,  fl  =  const.  Front  shock  maps  into 
/1(a),  and  left  and  right  sides  of  tail  shock  into  /lt(a)  and  /?3(<z),  respectively. 


where  A(a)  is  an  arbitrary  function  to  be  determined  later,  and  we  recall 


A  =  (y  +  i)/(y- 1). 

Similarly,  front  (5)  we  get 

»v  y(a,A)  =  y(M)+Jo  A(a) a~x sin  (0  +  /i) da. 


At  p  -  /1(a)  the  condition 


dx  yx  +  xp/l'(a) 


(15) 


(10) 


must  be  satisfied.  Elimination  of  y  using  (5)  and  substitution  of  (14)  for  x  produces  a 
linear  integral  equation  for  A(a): 


A(a.)Q(a, /}(<*)) +b(a)  1  +  A(a)Qp(a,/l(a))dd 


=  0, 


(1") 


where  Q  =  a~x  cos  (0  +  /i),  and 


6(a)  =  /?'(«) 


tan  r)  —  tan  0 
tan  ?/  —  tan  (0  +  ft) 


If  the  solution  for  0,  /i  and  s  is  known  in  the  a/1  plane,  this  equation  can  be  solved  for 
A(a),  and  the  transformation  back  to  the  physical  plane  computed  with  (14)  and  (15) 
In  general,  however,  the  solution  in  the  a/1  plane  depends  on  x,  through  (10). 

Up  to  this  point  the  equations  in  a(l  co-ordinates  have  been  derived  without  ap¬ 
proximation,  and  hence  are  equivalent  to  the  original  set  (2)  and  (3). 


V. 
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4.  Approximate  solutions 

If  the  third-order  changes  in  s  and  r~  are  neglected,  that  is,  it  is  assumed  that 
s  =  0  and  r~  =  -  P(p0)  everywhere,  the  solution  of  (2)  and  (3)  is  a  simple  wave,  in 
which  all  quantities  are  constant  on  the  principal  (C+)  characteristics,  which  in  turn 
are  straight  lines: 

0  =  tan-'/'(/?),  p  =  P~1(0+  P(p0)),  s  =  0  onC+: 

V  =/(/?)  +  (*-/?)  tan  (0+p). 

This  approximation  is  due  to  Friedrichs  (1948).  (Friedrichs  further  simplified  the 
problem  by  neglecting  terms  of  third  and  higher  order  throughout  the  calculation.) 

Because  simple  wave  theory  takes  s  and  r~  constant  at  their  upstream  values,  it 
can  be  expected  to  be  least  accurate  near  the  airfoil,  where  the  shock  is  strongest  and 
the  deviation  from  upstream  conditions  is  the  greatest.  An  improved  approximation 
in  this  region  can  be  obtained  using  shock  expansion  theory,  in  which  s  and  r~  are 
assumed  to  be  constant  at  their  values  just  behind  the  shock  at  the  leading  edge,  say 
s  =  s0  and  r~  =r„.  This  leads  to  a  slightly  modified  version  of  the  simple  wave 
solution : 

0  =  tan_1/(A).  /*  =  P~\0~ro),  s  =  s0. 

This  approximation  produces  a  very  accurate  solution  at  the  airfoil,  even  for  flows 
with  strong  shocks,  in  which  s  and  r~  are  not  at  all  constant  globally.  Hayes  &  Probstein 
(1966)  explain  that  the  down-running  waves,  which  can  be  considered  reflections  of 
the  outgoing  simple  wave  by  the  bow  shock,  are  fairly  weak  and  are  nearly  cancelled 
by  reflections  from  the  entropy  (or  vorticity)  layers.  Mahony  (1955)  gives  a  similar 
explanation.  The  shock  expansion  solution  rapidly  loses  accuracy  as  the  distance 
from  the  airfoil  increases.  This  is  in  contrast  to  simple  wave  theory,  which  is  more 
accurate  at  infinity. 

The  only  assumption  in  the  shock  expansion  solution  at  the  airfoil  is  that  r~  is 
constant.  Mahony  &  Skeat  (1955)  and  Meyer  (1957)  have  pointed  out  that,  since  any 
streamline  is  a  potential  airfoil,  r~  should  be  approximately  constant  along  each 
streamline,  that  is  r~  =  r~(a).  In  the  literature  this  assumption  has  been  employed 
in  various  ways.  If  r~  —  r~( a),  then  by  (10)  0  =  6(f),  i.e.  0  isconstanton  C+  characteri¬ 
stics.  This  in  turn  implies  that  the  pressure  is  constant  on  C+  characteristics,  as  can 
be  seen  from  the  following  form  of  (3  + ): 

d0  +  — -  —  =  0  onC'+:  ~  =  ta.n(d+p).  (18) 

2  y  p  dx 

Taking  both  0  =  0(f)  and  p  =  p(f)  along  with  r~  =  r~(a)  overdetermines  the  problem 
however,  since  any  one  of  0,  P  and  r~  can  be  written  as  a  function  of  the  other  two 
(and  s).  This  was  noted  by  Eggers  el  al.  (1953).  In  their  generalized  shock  expansion 
method  it  is  resolved  by  averaging  results  assuming  r~  =  r~(a)  and  0  =  0(f)  with 
those  assuming  r~  =  r~(a)  and  p  =  p(ft)  (see  Hayes  &  Probstein  1966,  p.  498).  Meyer 
(1957),  on  the  other  hand,  implicitly  drops  the  assumption  p  =  p(f),  and  uses  the 
solution  r  =  r'(a)  and  0  =  0(f),  which  satisfies  (10)  exactly,  but  does  not  satisfy  (7). 

In  the  present,  formulation,  it  appears  to  be  more  consistent  to  approach  the  prob¬ 
lem  in  either  of  two  ways:  in  equation  (10)  assume  either  (i)  the  left-hand  side  or  (ii) 
the  right  hand  side  is  zero.  Then  solve  (10)  together  with  the  remaining  equation,  (7). 
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In  case  (i),  the  solution  becomes  0  ~  p  =  p{fi)  and  s  =  s(a).  The  function  0{fi) 
is  determined  by  the  boundary  condition,  and  p(p)  must  be  determined  by  the  shock 
conditions.  It  then  happens  that  over  the  rear  half  of  the  airfoil,  /?  >  /?( 1 ),  p(fi)  cannot 
be  found,  since  no  data  is  specified  on  the  rear  shock.  This  difficulty  does  not  arise  in 
approach  (ii),  which  is  the  one  we  adopt. 

This  approach  can  be  thought  of  more  simply  as  arising  from  the  assumption  that 
0  is  constant  on  C +  characteristics,  rather  than  the  assumption  that  r~  is  constant  on 
streamlines.  If  0a  =  0,  then  (10)  reduces  to 

{0-P(/i))„  =  0  or  O-P(n)  =  -  /;(<*),  (19) 

where  P0(a)  =  P[ji(a, /?(a))]  -0(a,/?(a)).  Substitution  oft?  =  P(p)  -  P0(a.)  in  the  remain¬ 
ing  equation,  (7),  then  gives 

2P(,i)a-P0'(oc)  =  S-^s'(cc).  (20) 

P0(a)  and  s(a)  are  both  given  explicitly  by  the  shock  conditions,  so  (20)  can  be  re¬ 
garded  as  an  ordinary  differential  equation  for//,  in  which  /?  enters  only  as  a  parameter. 
It  is  nonlinear,  but  can  be  readily  solved  using  standard  numerical  methods.  The 
initial  and  final  values  of  /<  along  a  given  C r  characteristic  are  bot  h  given,  by  the 
boundary  condition  and  the  shock  conditions,  respectively,  which  allows  us  to  solve 
for  the  free  boundary  /?(a).  The  solution  in  the  a/?  plane  is  then  completed  by  com¬ 
puting  0(a,P)  =  P(ji(a,P))  —  P0(ol).  The  solution  for  0,  /t  and  s  in  the  a/?  plane  is 
independent  of  x  and  y,  because  (10),  the  only  equation  in  which  x  or  y  appears,  is 
neglected.  The  transformation  back  to  the  x y  plane  is  found  by  solving  (17)  for  A(a) 
(also  a  simple  numerical  calculation)  and  evaluating  the  integrals  (14)  and  (15).  The 
solution  obtained  from  this  approximation  will  satisfy  the  boundary  condition  and 
all  three  shock  conditions,  but  will  satisfy  (10)  only  approximately. 

This  approach  requires  more  work  (the  solution  of  an  ordinary  differential  equation 
on  each  C+  characteristic)  than  approach  (i)  or  the  generalized  shock  expansion  method, 
but  has  been  found  to  be  more  accurate.  Additional  support  for  this  choice  is  lent  by 
the  fact  that  the  factor  multiplying  0a  in  (10)  is  in  general  quite  small.  Approach  (i) 
has  however  been  found  useful  for  calculating  the  flow  behind  the  tail  shock,  where 
method  (ii)  is  difficult  to  employ  (see  appendix  B). 


5.  Numerical  method 

Our  approximate  solution  docs  not  satisfy  (Id),  or,  equivalently,  the  C~  equation 
(8).  In  this  section  we  present  a  simple  iterative  method  for  correcting  the  solution  so 
that  it  will  satisfy  all  the  equations  and  conditions. 

The  approximate  solution  is  computed  on  a  rectangular  grid  in  the  aft  plane  (as 
shown  in  figure  3),  which  is  then  vised  in  the  numerical  method.  The  front  shock  /9(a) 
is  therefore  kept  fixed  throughout  the  iterations.  This  fixes  (lie  normalization  of  a, 
so  for  every  iteration  beyond  the  original  approximation  J/(a)  is  not  given  by  (12) 
and  must  be  found  as  part  of  the  solution.  This  also  implies  that  a  =  1  no  longer  will 
correspond  exactly  to  x,  y  -»  oo. 

Given  the  approximate  solution  for  0,  //,  s  and  r  in  the  a/?  plane,  a  corrected  value 
of  r~  is  computed  from  the  C~  equation  (3—  ),  or  (S),  starting  at  the  shock  with  the 
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value  given  by  the  shock  conditions  and  numerically  integrating  downward  along  the 
C~  characteristics: 

_  f  8in2/t 

*8hOCk  2  y  (21) 

In  particular,  this  determines  a  new  value  r~(0,ft)  at  the  airfoil,  which  determines 
a  new  value  of  r4( 0,ft)  there,  since  r+  =  2 0  —  r~,  and  0(0,  ft)  is  given  by  the  boundary 
condition.  With  this  as  an  initial  value,  a  new  r+  is  computed  everywhere  by  numeri¬ 
cally  integrating  (3  +  ),  or  (7),  along  C+  characteristics: 

r+(a,y?)  =  r+(0,ft)+j‘?^s'(cc)d<x.  (22) 

With  r4  and  r~  thus  determined,  the  solution  given  by 

0=|(r+  +  r-),  p  =  P_I{|(r+ — r-)], 

and  s  will  satisfy  the  differential  equations  and  the  boundary  condition.  However,  the 
new  value  of  r4  (a, /7(a))  will  not  in  general  satisfy  the  shock  conditions,  and  hence 
will  imply  a  different  value  for  the  shock  angle  ij(a).  This  can  be  used  to  determine  a 
new  initial  value  r~(a,ft(a))  for  integrating  (21),  and  the  procedure  can  be  repeated. 

The  transformation  back  to  the  xy  plane  is  found  by  numerically  solving  the  integral 
equation  (17)  and  evaluating  the  integrals  (14)  and  (15).  This  must  be  done  at  each 
iteration,  since  x  and  y  enter  into  the  computation  of  the  integral  in  (21).  The  C~ 
characteristics  are  oblique  to  the  (a,/7)  co-ordinate  system,  so  at  each  point  a  small 
section  of  the  C~  characteristic  through  that  point  is  extended  backwards  to  intersect 
a  grid  line,  and  a  one-step  integration  is  used  to  compute  r~.  We  might,  in  place  of 
equation  (8),  have  integrated  (10),  which  has  the  advantage  that  r~  is  differentiated 
only  with  respect  to  ft,  so  that  the  integration  would  be  along  the  co-ordinate  lines, 
as  in  (22).  In  practice,  however,  this  has  been  found  unadvantageous.  The  solution  does 
not  converge  as  quickly,  and  may  not  converge  at  all  without  modification  (see  Chong 
&  Sirovich  1980).  We  attribute  this  to  the  fact  that  small  variations  in  r~  are  naturally 
propagated  along  the  C~  characteristics. 

This  scheme  has  been  implemented  using  second-order  numerical  methods  (trape¬ 
zoidal  rule,  improved  Euler  method,  etc.).  Some  results  are  given  in  the  next  section. 


6.  Results 

Calculations  have  been  performed  for  several  airfoils  over  a  range  of  Mach  numbers. 
The  results  presented  in  figures  1  and  3-7  are  for  a  symmetric  circular  arc  airfoil  with 
thickness  ratio  0-25  at  upstream  Mach  number  Ma  =  4.  In  figures  8-10  results  from 
the  additional  cases  M0  =  2-5  and  7-5,  for  the  same  airfoil,  are  included  sis  well.  These 
cases  were  chosen  in  part  for  the  interesting  effects  they  exhibit. 

The  iteration  scheme  converges  quite  rapidly,  based  on  a  comparison  of  the  solutions 
at  successive  iterations.  In  table  1,  the  maxima  (over  all  grid  points)  of  the  differences 
in  the  values  of  0,  /i  and  x  are  given  for  the  case  M0  =  7-5  (the  most  slowly  convergent 
of  the  three  cases).  The  greatest  differences  are  in  x  and  usually  occur  near  a  =  1, 
where  x  ->  oo.  The  errors  in  x  are  smaller  closer  to  the  airfoil.  For  thinner  airfoils  or 
lower  Mach  numbers,  fewer  iterations  arc  required  for  the  same  accuracy.  In  the  case 


Iteration 

A0/0(0,  0) 

A/i/p 

A.r/x 

1 

o-05»r> 

0  1 170 

0  3701 

2 

0  0141 

0  0000 

0  1221 

3 

0  0008 

00007 

0-0112 

4 

0  0002 

0  0000 

0  0017 

5 

00001 

0  0003 

0  0000 

of  a  10%  thick  parabolic  arc  airfoil,  for  example,  even  at  M0  =  10  the  difference 
between  the  approximate  and  exact  solutions  is  less  than  one  per  cent  in  0  and  p  and 
six  per  cent  in  a\  In  such  a  case  there  is  little  reason  to  go  beyond  the  approximate 
solution. 

The  case  M0  =  2-5  is  discussed  in  Holt  ( 1 077).  Figure  4  contains  a  comparison  of  the 
leading  shock  when  computed  by  our  approximate  and  exact  methods,  the  BVLR 
method  (an  exact  numerical  method),  and  generalized  shock  expansion  theory  (the 
latter  and  the  BVLR  solution  are  taken  from  Holt  1077,  p.  77).  In  this  case,  our 
approximate  solution  is  indistinguishable  from  the  exact  solution.  The  small  difference 
between  these  and  the  BVLIl  solution  is  probably  attributable  to  copying  errors. 

Figure  5  contains  plots  of  pressure  contours  in  the  xy  plane  and  the  value  of  log p 
on  the  airfoil  surface  and  on  the  line  of  symmetry  behind  the  airfoil.  Comparison  with 
figure  1  shows  that  the  contour  lines  between  the  lead  and  tail  shocks  are  nearly 
identical  to  C+  characteristics,  i.e.  the  pressure  is  approximately  constant  on  C ~ 
characteristics.  This  was  seen  in  §4  to  be  related  to  the  fact  that  0  is  approximately 
constant  on  C+  characteristics,  which  in  turn  is  related  to  the  fact  that  r~  is  approxi¬ 
mately  constant  on  streamlines.  The  latter  two  assumptions  are  illustrated  in  figures 
6  and  7. 

In  figure  6,  the  deflection  angle  0  is  plotted  versus  a  on  each  of  the  C+  characteristics 
shown  in  figure  3.  In  the  region  behind  the  tail  shock  0  is  very  nearly  zero  (|d|  <  0-005) 
everywhere.  The  variation  in  0  along  each  characteristic  is  quite  small,  with  the  most 
serious  departure  occurring  on  the  characteristics  originating  from  the  rear  part  of 
the  airfoil.  These  characteristics  tend  to  intersect  the  tail  shock  fairly  close  to  the  air¬ 
foil,  however.  A  related  phenomenon  is  that  the  principal  characteristics  are  nearly 
straight.  This  however  does  not  remain  true  in  the  region  behind  the  airfoil. 

Figure  7  shows  the  variation  of  r~  with  /?  on  each  streamline  of  figure  3.  Somewhat 
remarkably  the  assumption  r~  =  —  P0(a)  is  bet  ter  at  the  airfoil  than  a  short  distance 
away.  The  assumption  is  less  satisfactory  behind  the  tail  shock.  The  rapid  downstroke 
of  the  r~  curves  also  indicates  a  large  value  of  0X,  although  0  itself  remains  quite  small. 

The  entropy  jumps  created  by  the  lead  and  tail  shocks  are  given  in  figure  8  for  the 
three  cases  M0  =  2-5,  4-0  and  7-5.  The  entropy  variation  along  the  tail  shock  has  a 
two-scale  appearance,  especially  at  the  higher  Mach  numbers,  which  shows  a  very 
rapid  decrease  in  strength  in  the  initial  portion  of  the  shock.  The  slower  variation  in 
entropy  follows  that  induced  by  the.  front  shock.  looking  at  figure  1,  we  see  that  the 
streamlines  spread  apart  rapidly  as  the  flow  passes  the  midchord  position.  The  incli¬ 
nation  of  the  flow  incident  upon  the  tail  shock  therefore  decreases  rapidly,  which 
causes  a  correspondingly  rapid  decrease  in  shock  strength. 

Another  important  effect  is  also  at  work  in  this  region.  The  gas,  which  is  compressed 
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Ficukk  4.  Front  shock  for  flow  field  of  figure  1  as  computed  by :  present  approximate  and  exact 
methods  ( - ),  BVLR  method  ( - ),  generalized  shock  expansion  method  ( - ). 

at  the  front  shock,  in  following  the  profile  past  the  midchord  experiences  a  rapid 
expansion,  which  is  strong  enough  that  the  local  Mach  number  at  the  trailing  edge 
exceeds  the  upstream  value  (M  =  9-83  for  the  M0  =  7-5  case).  This  recovery  process 
is  largely  cut  off  by  the  tail  shock,  however,  since  the  large  negative  value  of  0  on  the 
after  part  of  the  airfoil  causes  the  principal  characteristics  to  have  negative  slopes, 
so  that  waves  originating  there  must  intersect  the  tail  shock  near  the  airfoil.  As  a 
result  the  Mach  number  along  the  tail  shock  falls  off  rapidly,  which  augments  the 
rapid  decrease  in  strength  of  the  tail  shock.  For  the  case  Jl/0  =  7-5  the  Mach  number 
along  the  shock  even  falls  below  7-5. 

The  pressure  field  behind  the  airfoil  (figure  5)  also  contains  interesting  features.  In 
spite  of  the  very  high  shock  strength  at  the  trailing  edge,  the  pressure  jump  through 
the  shock  does  not  quite  bring  p  up  to  the  equilibrium  pressure  p  =  1 .  There  is  a  rapid 
pressure  increase  immediately  behind  the  trailing  edge,  in  which  p  increases  above  the 
equilibrium  value,  reaching  a  maximum  about  one  chord  length  out.  The  return  to 
equilibrium  from  this  point  is  very  gradual.  The  total  variation  in  pressure  behind  the 
tail  shock  is  quite  small  compared  with  that  along  the  airfoil  surfaces. 

Far  behind  the  airfoil  p  -*■  1  and  0  -*■  0.  It  then  follows  from  the  equation  of  state 
that 

n2  =  exp  ( —  (y  —  1 )  53(a)/ y  ], 

where  s3(a)  is  given  by  figure  8.  From  (1),  we  can  then  compute  the  velocity  q  at 
infinity.  This  is  shown  in  figure  9  for  M0  =  2-5,  4  0  and  7-5.  As  a  result  of  the  non- 
uniform  entropy,  the  flow  at  infinity  has  a  vorticity  distribution. 

A  feature  which  is  difficult  to  perceive  from  figure  1  or  figure  5  is  that  the  tail  shock 
angle  is  not  monotonic.  In  figure  10  the  variation  of  the  sloj)C  of  the  tail  shock  is 
given  for  the  three  cases  we  have  discussed.  In  each  case  the  shock  angle  decreases 
on  leaving  the  trailing  edge.  (This  result  has  been  verified  indej>endcntly  by  J.  ('. 
Townsend  1979  (private  communication),  using  a  numerical  method  developed  by 
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Figure  5.  Upper  graph:  pressure  contours  in  flow  field  of  figure  1.  Lower  graph  : 
log  pvs.x  at  a  —  0.  Note  different  scales  for  0  <  .r  <  1  and  .r  >  1 . 

M.  D.  Salas.)  This  is  contrary  to  what  is  observed  for  lower  Mach  numbers  or  thinner 
bodies.  We  have  seen  that  the  inclination  of  the  incident  flow  decreases  along  the 
shock.  If  the  Mach  number  upstream  of  the  shock  were  constant,  this  would  predict 
a  decrease  in  shock  angle.  The  Mach  number  actually  decreases  along  the  shock  how¬ 
ever,  which  tends  to  increase  the  shock  angle.  At  high  Mach  numbers  the  shock 
angle  is  more  dependent  on  the  flow  angle  than  on  the  Mach  number,  as  can  be  seen 
from  the  fact  that  the  shock  polars  for  different  Mach  numbers  approach  a  limiting 
curve  as  M  -*■  oo  (see  e.g.  Liepmann  &  Roshko  1957,  p.  87).  In  these  cases,  near  the 
trailing  edge  the  decreasing  flow  angle  dominates.  Farther  away  from  the  airfoil,  or 
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Figure  10.  Tail  shock  slope  tun  i/,(a)  for  M0  -  2  f>,  4  0  and  7 
Dashed  lines  are  asymptotic  values,  tan  //„. 


in  problems  with  lower  Mach  numbers  or  thinner  airfoils,  the  effect  of  decreasing  Mach 
number  dominates. 

In  the  case  jl/0  =  1-5  the  shock  angle  undergoes  a  second  oscillation  in  which  it  rises 
above  the  Mach  angle  at  infinity,  /i0.  This  is  explained  by  the  rapid  fa"-ofT of  Mach 
number  along  the  shock,  below  its  value  at  infinity.  A  final  item  of  note  in  figure  10 
is  that  for  M0  =  7-5  the  shock  angle  actually  starts  off  with  a  value  which  is  greater 
than  /<„.  As  il/0  ->  oo  the  upstream  Mach  angle  //„  goes  to  zero,  as  does  the  Mach  angle 
at  the  trailing  edge,  since  the  Mach  number  there  also  increases.  The  shock  slope  at 
the  trailing  edge  approaches  a  finite  value  however,  which  de|>cnds  on  the  airfoil 
slope  at  the  trailing  edge. 


7.  Conclusions 

The  methods  we  have  presented  are  useful  in  computing  two  dimensional  flow  fields 
about  airfoils.  The  approximate  solut  ion  is  accurate  enough  for  many  cases  of  interest . 
and  the  numerical  method  furnishes  a  rapid  correction  to  the  solution  in  those  cases 
where  it  is  not.  The  characteristic-streamline  co-ordinate  system  is  useful  both  >r 
the  computation  of  the  approximate  solution  anti  the  corrections,  and  is  also  e 
venient  for  displaying  and  interpreting  the  results. 
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The  use  of  (he  streamlines  as  one  co-ordinate  and  the  iterative  nature  of  the  numeri¬ 
cal  calculation  make  the  met  hod  convenient  for  the  incorporation  of  a  boundary-layer 
correction.  In  a  boundary  thickness  method,  for  example,  a  succession  of  inviscid 
calculations  are  performed  with  a  changing  airfoil  shape.  The  changing  shape  could 
he  easily  included  in  the  present  iteration  method. 

In  response,  to  a  referee’s  request  for  comparison  with  other  integration  schemes, 
we  asked  I)r  dames  C.  Townsend  of  the,  NASA  Langley  Research  Center  to  run  some 
speed  trials  on  their  01)0  Cyber  175  computer  comparing  our  code  with  a  ‘marching' 
method  developed  there.  At  the  lowest  Mach  number,  Mn  =  1-25,  our  scheme  runs 
about  seven  times  faster  than  the  marching  method,  while  at  the  highest  Mach 
number,  Ma  =  10,  our  scheme  was  slightly  slower.  The  present  method  is  most  efficient 
at  low  Mach  numbers  where  the  approximate  solution  is  most  accurate  and  the  fewest 
iterations  arc  required.  This  is  in  contrast  to  the  marching  method,  where  low  Mach 
number  necessitates  a  short  step  size  for  stability,  and  hence  longer  computation  times. 
While  these  trials  give  some  idea  of  relative  speed  they  cannot  be  considered  definitive. 

This  work  was  supported  by  the  National  Aeronautics  and  Space  Administration 
under  NASA  Grant  no.  NSG  l(il7.  The  authors  would  like  to  thank  Dr  James  0. 
Townsend  for  carrying  out  a  number  of  computations  which  were  very  useful  in  the 
course  of  this  research. 

Appendix  A.  Case  of  an  arbitrary  gas 

For  an  arbit  rary  gas,  the  equations  of  motion  in  characteristic  form  can  be  written 
(Haves  &  Prohstein  19G(>,  p.  484) 

ds  —  0  on  dy/dx  =  tariff,  (A  1) 

d0±  <1  >dp  =  0  on  dy/dx  =  tan  (0  +  //),  (A  2) 

where  <1>  =  pj(i>aa\pq~  tan/r).  We  can  consider  O  to  be  a  function  of  p  and  s.  By 
introducing  the  variables 

<o(p,  .$)  =  J  <l>(p,  s)dp  and  Q(p,s)  =  dcj(p,s)/cJ$, 
which  are  defined  so  that  dw  =  <l>dp  +  Qds,  (A  2)  can  be  written  as 

dO±dio  =  +  Lids  on  dy/dx  =  tan  (ff  ±/i).  (A  3) 

If  o)  and  12  are  now  regarded  as  functions  of //  and  s,  (A  1 )  and  (A  3)  are  three  equations 
....  .  in  three  unknowns:  ff,  p  and  s.  Equations  (3)  arc  a  special  case  of  (A  3)  in  which 

m  =  /»(//)  and  12  =  (sin  2//)/2y. 

The  transformation  to  a/1  co-ordinates  goes  through  for  the  most  part  as  before. 
Equations  (<>)  (S)  in  the  general  case  become 

•S'/i  ;  d,  (fl  +  Oi):  -=  12.s '(a), 

where  ir  is  si  ill  given  bv  (!i).  The  counterpart  of  ( 1 3)  is 

<»  -t  (/i  t  fuj^cot//  f  (ff  +  /<)/i tan  (0  +  /i ). 
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This  equation  can  in  principle  be  solved  in  the  same  manner  as  (13),  but,  depending 
on  the  form  of  we  may  not  have  an  explicit  integral  like  (14). 

The  assumption  0a  =  0  in  the  general  case  implies  (0  —  =  0  or  8  —  w  =  —  w0(a). 

The  resulting  approximation  can  be  expected  to  be  valid  at  least  in  cases  in  which  the 
behaviour  of  the  gas  does  not  differ  too  greatly  from  that  of  a  perfect  gas  with  constant 
specific  heats  and  y  =  1-4.  In  particular,  it  has  been  shown  (see  Hayes  &  Probstein 
I960,  §7.2)  that  shock  expansion  theory  tends  to  lose  accuracy  if  y  is  allowed  to 
approach  1. 


Appendix  B.  Tail  shock  for  a  symmetric  airfoil 

In  general,  the  solutions  above  and  below  the  airfoil  can  be  computed  independently, 
up  to  the  appearance  of  the  tail  shocks.  The  flows  from  the  top  and  bottom  interact 
behind  the  airfoil,  which  complicates  the  computation  of  the  tail  shocks  and  the  flow 
behind  them.  The  upper  and  lower  regions  behind  the  airfoil  are  separated  by  a 
contact  discontinuity,  or  slipstream,  across  which  0  and  p  are  continuous,  but  the 
other  variables  jump.  In  the  case  of  an  airfoil  symmetric  with  respect  to  the  x  axis 
the  slipstream  coincides  with  the  x  axis,  and  can  be  considered  a  rigid  boundary. 
The  problem  is  still  quite  different  from  the  front  shock  problem,  because  the  flow 
upstream  of  the  tail  shock  is  not  uniform. 

The  transformation  to  a/?  co-ordinates  behind  the  tail  shock  can  be  chosen  differently 
than  that  ahead  of  it.  In  particular,  it  is  more  proper  to  regard  the  C~  characteristics 
as  the  principal  characteristics,  since  the  C+  waves  are  only  produced  as  reflections 
of  the  C~  waves,  which  originate  at  the  tail  shock.  The  approximate  solution  is  some¬ 
what  more  accurate  if  the  C~  characteristics  are  used.  On  the  other  hand,  for  numerical 
work  it  is  better  to  take  the  C+  characteristics  as  the  /?  co-ordinates,  because  this  has 
the  effect  of  putting  more  points  near  the  trailing  edge,  where  a  rapid  variation  in  the 
solution  occurs.  We  keep  a  constant  on  streamlines  as  they  cross  the  shock,  and 
normalize  (i  behind  the  tail  shock  so  that  the  infinite  region  behind  the  tail  shock  is 
mapped  into  a  finite  region  in  the  a/l  plane.  In  the  calculations  presented  here,  this 
was  done  by  setting  /?3(a)  =  1  +  £a,  producing  the  triangular  region  shown  in  figure  3. 

The  approximate  solution  used  for  the  flow  over  the  airfoil  cannot  be  conveniently 
employed  for  the  flow  behind  the  tail  shock,  because  the  non-uniform  flow  to  its  left 
makes  it  impossible  to  calculate  /’0(a)  and  s(a)  a  priori  for  use  in  (20).  Therefore  the 
simpler  of  the  approximations  given  in  §  4  is  used:  0  =  03(P),  p  =  p3(p),  and  s  =  s3(a). 
All  the  characteristics  intersect  the  x  axis,  where  0  =  0,  so  03(P)  =  0,  and  hence  in  this 
approximation  0  =  0  everywhere.  This  turns  out  to  be  quite  accurate  (see  §6).  Given 
that  0  =  0  behind  the  tail  shock,  it  is  possible  to  solve  the  shock  conditions  for  the  tail 
shock  angle  rj2(a),  in  terms  of  the  solution  upstream  of  the  tail  shock,  which  we 
assume  has  been  previously  computed.  This  also  determines  p3(P)  and  s3(a),  and  gives 
an  ordinary  differential  equation  to  solve  for  the  tail  shock  /?2(a).  It  is  possible  to 
derive  expressions  for  x  and  y  similar  to  (14)  and  (15)  for  the  region  behind  the  tail 
shock,  which  will  involve  a  new  function  A3(a).  An  explicit  solution  for  A3(a)  can  be 
found  in  this  case,  involving  the  computed  tail  shock  trajectory. 

The  iteration  scheme  proceeds  essentially  as  before.  Given  r~(a,/?3(a))  from  the 
shock  conditions,  we  integrate  (21)  along  C~  characteristics  down  to  the  slipstream 
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a  =  0.  Then  we  reset  r*( 0,fi)  =  -r~(0,/}),  and  integrate  (22)  upwards  to  /?3(a).  The 
new  r+  and  r~  define  a  new  0(a.,/i3(ct)),  which  is  used  to  solve  for  a  new  shock  /?3(a)  and 
new  functions  ?/2(a),  s3(a),  and  r~(a,/3s(a)),  with  which  we  start  the  next  iteration. 
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Technical  Notes _ 
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Introduction 

III-  inverse  oi  design  problem  lot  the  case  ol  iwo- 
dimension. il  stipci  sonic  an  toil  shapes  is  considered.  In 
view  of  the  hypcihohc  stinelme  ol  die  imdciking  evpiaiions. 
die  calculation  is  simplei  than  the  corresponding  subsonic 
problem.  In  a  recent  paper.1  die  authors  developed  a 
numerical  procedure  loi  Heating  the  direct  problem,  the 
calculation  of  supersonic  flow  fields  past  erven  profiles.  This 
procedure  makes  use  ol  streamlines  as  one  ol  the  coordinates. 
\s  a  result,  il  is  especially  suited  to  the  inverse  problem.  Since 
the  adaptation  ol  the  method  to  the  present  problem  is  very 
similar  to  the  oriental  lormulalion,  vve  will  give  only  a  biief 
outline  ol  the  pioeediue  in  this  Note,  lor  purposes  of 
comparison,  vve  also  ptesenl  two  approximate  solutions  of  the 
problem,  one  a  simple  trealmenl  based  on  lineari/ed  analysis 
and  lire  other  based  on  shock  expansion  theory.  The  latter 
proves  to  be  highly  accurate  and  fails  only  in  extreme  eases. 

Outline  of  the  Method 

Consider  two-dimensional  supersonic  flow,  which  we 
describe  by  die  flow  deflection  ancle  »,  the  Mach  angle 
p  sin  '(|r.W).  and  v.  the  entropy  divided  by  the  gas  con¬ 
stant  R. 

I  lie  equations  ol  motion  in  characteristic  form  are" 
siu2ii  di 

dr  ’  r  els  on  (  "  :  -a  tan  (tt  t  p )  ( I ) 

2-,  dv 


dv  -  0  on  stieanilines:  tan"  (2) 

el  v 

and  where  /’  "t/’(p)  anel  /'(e).  the  I’randtl  ancle,  is 

ile lined  bv 

/’( 1 1 )  X  tan  '  ( .X  l.iit/i )  X  -  (•)  r  /)/(-,  / )  (A) 

We  introduce  a  cootdmaie  system  made  up  of  the  streamlines 
(n  const)  and  the  f"  eliaiaelet istics  (,(  const)  Hie 
er|iialions  then  become 

v  0  (4) 

si  n2/i 

r\  ,  v  (') 

-  T 
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slll2,i  2  \ 

r  +  nr  =  —  v  .  n  -  -  (ft) 

2',  I  -  tanUtanp  v 

In  this  fotimilation.  the  eooulinates  (v.e)  in  the  physical 
plane  become  dependent  variables  and  are  cover  lied  by  the 
equations 

i  =  v  t  an//.  r  ■ ,  -  v  tan(d  *  p )  f) 

liquations  (4-7)  are  to  be  augmented  by  the  shock  relations, 
which  are  not  repeated  here,  and  by  the  given  airfoil  presstne 
distribution  />  =  /)„(  v) . 

The  transformation  from  the  physical  to  the  (u.d)  plane 
leaves  open  two  arbitral y  Itrnelions.  One  ol  these  is  lived  by 
set  tint! 

\tt)..i)--.i  (SI 

on  the  streamline  <»  0.  which  wo  take  to  be  the  as  vei 

unknown  airfoil.  As  a  second  condition,  we  requite  dial  the 
shock  ancle  rj  vary  linearly  with  <>, 

i| (n  )  =  i) (/))  +  I  p„  •  <i(0)  | if.  0<  it  -  I  (9) 

The  method  of  solution  is  itetative  and  begins  with  an 
approximate  solution.  T  his  approximation  is  allied  to  shock 
expansion  theory.'  1  f  igure  I  contains  a  sketch  ol  the  plane  ot 
integration.  The  curve  n(d)  represents  the  as  yet  unknown 
shock  trajectory.  As  indicated  in  I  ig.  I,  a  uniform  n  mesh  n 
chosen,  which  with  nr  (if)  generates  the  d  mesh  over  the  lirst 
portion  of  the  figure.  The  determination  of  the  solution  starts 
with  the  replacement  of  lq.  (0)  by  the  approximation 

/•  =«-/>(,,)  =  -/>„(.,)  (1(1) 

where  -/*„(<«)  is  the  value  determined  at  the  shock  II  lq 
( 10)  is  substituted  into  I- q.  (5).  we  obtain 

0  sin2u 

2.  /'(/<  >  -/».:(<.)  +■  v(m)  (II) 

(In  2', 

w  hete  v  ( u  )  is  also  determined  from  the  shock  telat ions. 

lot  each  given  value  ol  .i.  this  is  an  ot  dinars  dilleienii.il 
equation  lot  p.  which  can  be  solved  numerically.  Initial  and 
final  values  p|.i(d),dl  and  u(0,d)  are  both  known,  which 
allows  us  to  determine  the  shock  n  (o' 1 .  1  he  values  p  |  o  I  >).  i 1 
come  from  the  shock  conditions  and  the  values  p(0, d)  Itom 
the  billowing  (elation  between  p.v.  and  the  not  m.ih/cd 
pressure  /' 

esp I  '  ( v  r  (.,/>)  I 

..■>-/  1  7  1 

sin  p  (  I  _  l 

’  V  /  .  I  -  /  I 

/  t  y  A  1-„  -  esp  ( v  ^  i„/i) 

Since  v(0)  is  known  and  /i(O.if)  />.(.>)  is  given,  iliis 
detei mines  p(0.d ) . 

On  the  portion  ol  >■  0  not  lying  undei  the  Itonl  shock,  we 

choose  a  uniform  o'  mesh  We  use  lq  (12)  to  deieiuune 
p((),,f)  and  integrate  lq  (II)  up  Irom  o  0  to  linislt  lire 
determination  ol  p(cr.,()  1  quation  (10)  is  used  to  deteimme 
0(  t. o').  I  lie  appi osimate  solution  is  then  known  evetywltete 
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in  i he  (nr.  d)  plane.  biiially.  bqx.  (7)  are  integrated  lo  t inti  the 
transformation  lo  die  physical  plane  everywhere.  In  par 
lieular.  die  hotly  shape /(.v)  is  given  by 


/<  v)  =/«))  +  |  ian»(«.i<)tlrf 


which  completes  die  eonipmaiion  of  die  approximate 
solution. 

I  lie  ilcralion  procedure  siarls  with  die  neglected  hqs.  (6), 
which  were  replaced  in  die  approximate  solmion  by  bq.  (10). 
We  numerically  integrate  b.qs.  (6)  downstream  along  the  (' 
eliaraeteristies,  starling  at  the  shock.  This  produces  new 
values  of  r  (it.ii)  everywhere.  Irom  /  (Ojf)  we  can 
determine 

r‘(0ji)=r  (O.ii)  +  2/’|, ■((/., f)  |  (14) 

since  /<(d.o')  is  given  by  bq.  (12).  liquation  (5)  is  then  in¬ 
tegrated  up  from  the  body  along  the  C '  eliaraeteristies  to  give 
/•  ■  (o.d)  everywhere.  We  can  then  obtain  improved  values  of 
(/and  1 1  everywhere  from 


1  I  Zi(r‘  -r  )| 


A  new  shock  angle  i)(nr)  is  determined  from  the  shock 
relations,  and  a  new  transformation  to  the  physical  plane  is 


obtained  from  I  qs.  ( I  J)  and  (7).  Hie  i  let  at  ion  is  then  repealed 
until  a  convergence ci itcrion  is  met. 

Approximate  Calculations 

I  inijri/iil  Approximation 

Perhaps  the  simplest  ealeulation  for  /(  v)  given  the  picssutc 
distribution  /i„(.v)  follows  from  linearized  llieot  v . '  In  out 
normalization,  lincaiizcd  theory 

s  Mi,  -  / 

/'  (  v)  =  l/>„(  v  )  -  /  | 

7  V  / , ; 


S  ,\fi  -  /  /  f  '  l 

/<  x >  =/(«)+  (\  />„(ii di-.v) 

"t  '•n  •  " 


shock  Kxpansion  Approximation 

One  max  also  base  an  approximate  ealeulation  ol  /( v )  on 
shock  expansion  theory.1  We  recall  that  if  (  i )  and  /■  (  v  ) 
denote  values  on  the  airfoil,  then  in  tilts  approximation  it  is 
assumed  that 

(/„  (  v)  -  l’\ ii„  (  v )  |  /'„  (17) 

where  is  the  value  at  the  leading  edge  behind  the  shock 
Lquation  (12)  detet mines  u„(v)  Irom  the  given  picssutc 
/> ,  (  v)  and  the  entropy  at  the  leading  edge  v„,  .mil  then  /(  v  l 
tollovvs  from 


/  '  - 


/;  v)  -  /((/)  i  |  t.in|/’( /.„  ( i. )  )  t’J'ii 


•  \  i 

!  : 

/f  :  i  :  ■ 

/..  :  .L.  i... 


■  •  l_:. 


!  i  i 

i  :  .. 


: ;  •  •  rrt|- 


As  an  examination  shmu,  /(  \  )  calculated  in  ihis  \\a\  apices 
with  the  li  r  si  approximation  in  the  iterative  volution  out  lined 
m  the  pre\  ions  section. 


I  i^.  I  Plane  of  integration:  lines  a  =  const  represent  streamlines, 
lines  ,i  const  represent  C  ’  characteristics,  anil  curve  represents 
the  l»o«  vhock. 


I  in.  A  Airfoil  shapes  computed  lor  M„  2.5  and  airtoit  pressure 
dislrilmtion  !.«/», ,(  v  1  I  2a. 


l-itt.  2  Airtoit  stripes  computed  for  -  1.2  and  airfoil  pressure 
distribution  0.05(1  -  2jtI * . 


I  ij».  4  Airfoil  shapes  eompuled  lor  \f„  5  and  airfoil  pressure 

distribution  i. >/), ,( t )  5(1  2\i. 
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Results 

Sample  calculations  are  shown  in  Figs.  2-4.  In  each  the 
results  obtained  from  linearized  theory,  shock  expansion 
theory,  and  the  exact  numerical  calculation  are  compared. 

Figure  2  contains  the  results  for  a  low  Mach  number 
(M0  =  1.2)  and  small  pressure  jump.  As  should  be  expected, 
all  of  the  results  are  in  close  agreement.  For  Fig.  3,  the  Mach 
number  is  moderate  (M„  =  2.5)  and  the  jump  in  t*p  at  the 
leading  edge  is  unity.  In  this  case,  linear  theory  is  poor  in 
predicting  an  overly  thick  body.  The  result  based  on  shock 
expansion  theory,  on  the  other  hand,  is  virtually  in¬ 
distinguishable  from  the  exact  case.  In  the  final  example  (Fig. 
4),  both  the  upstream  Mach  number  ( M0  =  5)  and  the  pressure 
jump  are  relatively  large.  Linearized  theory  is  now  very  poor. 
Shock  expansion  theory  still  does  quite  well  for  most  of  the 
derived  airfoil  and  begins  to  depart  significantly  only  near  the 
trailing  edge. 

Conclusions 

A  method  for  the  design  of  two-dimensional  supersonic 
airfoils  has  been  presented,  which  incorporates  available 
physical  and  mathematical  knowledge  of  the  problem  (e  g., 
shock  expansion  theory  and  characteristics),  in  order  to 
facilitate  the  numerical  computation.  A  similar  approach 
should  prove  useful  in  the  more  complicated  problem  of  the 
design  of  real  airfoils  in  which  three-dimensional  flow, 
boundary-layer  effects,  etc.,  rust  be  considered.  The  iterative 
nature  of  the  present  method,  in  particular,  makes  it  well 
suited  to  the  inclusion  of  boundary-layer  corrections. 
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by  1 1,  the  position  of  the  shock  is  governed  by 
dr 

tann  (3) 

dv 
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Brown  University, 
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We  introduce  new  coordinates  (ir.d)  ihrough 

r  v  lanM  (4) 


Introduction 

SUPl  RSONIC  inviscid  flow  can  generally  be  solved  by  the 
method  of  characteristics  or  by  shockcapturing 
methods,1  '  The  method  ol  characteristics  computes  the  floss 
along  characteristics  and  uses  the  Rankine-Hugoniot  condi 
lions  at  the  shock.  This  method  has  the  advantage  ol  ac¬ 
curacy,  but  is  regarded  as  complex  and  computationally  incl 
ficient,  especially  in  regions  of  near  coalescence  of  the  two 
sets  of  characteristics11  In  shock  capturing  methods  the 
shock  is  smeared  over  several  grid  points,  where  oscillations 
can  occur  and  the  scheme  loses  accuracy  However,  due  to 
their  directness  and  computational  ease,  shock  capturing 
methods  have  been  preferred  in  recent  seats 

In  this  Note,  we  develop  an  efficient  method  using  the 
characteristics  and  streamlines  of  a  flow  I  hose  are  used  as 
coordinates  and  the  floss  quantities  are  expressed  m  terms  ol 
Riemat.n  functions.  A  scheme  is  obtained  which  is  signtfi 
candy  more  efficient  and  accurate  than  shock  capturing 
methods  for  How  over  axisymmetric  bodies  Since  streamlines 
form  one  of  the  coordinates,  we  naturally  obtain  a  body-fit 
system.  It  is  also  a  truly  shock-fit  coordinate  system  Not  only 
are  the  Rankine-Hugoniot  conditions  used,  hut  the  shock  lies 
exactly  on  grid  points  also. 

rite  success  of  the  present  method  in  (he  two-dimensional 
case  rests  on  the  discovery  of  an  accurate  and  simple  approx 
until  ion.*  In  the  axisymmetric  case,  a  similaily  accurate  and 
simple  approximation  has  eluded  us  Hie  approximation 
presented  herein  is  simple,  hut  generally  not  as  accuiate  as 
that  for  the  two-dimensional  case.  A  better  iterative  pto 
cedure  has  been  developed  to  compensate  lor  (his  weakness. 
As  was  the  case  lor  the  two-dimensional  flow,  our  method 
is  well  suited  for  the  inverse  design  problem  (i.e  .  given 
AT,,  >  I  and  the  pressure  distribution  on  the  body,  find  the 
shape  of  the  body  and  the  flow  everywhere) 

Formulation 

VV  e  consider  axisymmetric  flow  with  incident  Mach 
number  A/,,  I  and  shock  attached  at  the  tip  I  he 
chaiaclenslic  equations  can  be  written  m  (cults  ol  entropy  s. 
Moss  angle  H.  Mach  angle  p  sin  1  (I  A/l  t  W  local  Mach 
number),  pressure  />,  density  and  velocity  q  as  follows* 

dr 

dv  0  on  tariff  ( I ) 

dv 

slllltsl  ll/i  dr  dr 

t  on  (  ’  (auto r  1 1 ) 

sin(M  i  /, )  r  d  v 

(3) 

At  the  hods  r  / (  x )  we  impose  the  boundaix  condition 
I.iiiM  /  (\)  lumps  acioss  a  shock  are  given  hx  the 
Rankine  Hugoniol  conditions  "  If  we  denote  the  shock  angle 


dn  i 


d/' 

in/  tnii/i 
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fig  1  Rolls  and  shuck  nt  Atlfl o-lhich  parahntu  Imrlv  ai  M„  2 


fig.  2  Pressure  distributions  on  the  body  for  I0*’i  and  AO*r»- thick 
parahoti,  bodies  al  Mn  2. 
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Here  y  -  1.4.  An  allernale  form  of  lq.  (7)  used  laier  is  as 
follows: 

K,=(l  ■  lanflian/i)-'-’-#  +  lan/i  (8) 

.v„  r 

ll  is  imporlanl  to  note  that  the  dynamical  equations  must  he 
augmented  by  the  coordinate  equations  (4)  and  (5)  I  Ins 
transformation  is  still  undetermined  up  to  two  arbitrate 
functions.  We  fix  one  of  these  functions  by  taking  v(0,/()  =  /f 
at  lhe  body  nr-O.  It  therefore  follows  dial 


dimensional  and,  hence,  shock  expansion  theory  ( K  con 
stain  along  an  entire  si  cam  I  me)  is  valid  Indeed,  tor  hyper 
some  llow,  along  f"  we  have  s  (a)  huge  and  r  small  so 
dial  Iqs.  (7)  and  (8)  aie  well  approximated  by  two 
dimensional  theory.  Al  lowei  Mach  numbers,  this  approx 
imation  is  no  longer  valid  However,  in  our  matching 
scheme,  we  need  only  to  approximate  one  grid  point  away, 
and  to  assume  that  R  changes  slowly  along  a  streamline. 
Combining  I  iqs.  (7)  and  (9). 


r(0.d)=/(d>.  fMO.fi)  =  tan  1  /  ( ) 


R  v  (I  tanfftaii/i) 


sin2|i  (J , 


To  determine  the  second  arbitrary  function  we  fix  the  shock 

b>  do 

.  .1  (101 


'  ,„(li  •  '  ) 

y  l  V  2xin  ,,  / 


l  he  above  loruiulaiion  must  be  modilied  slightly  in  order 
to  neat  the  inverse  design  problem  Tor  dux  problem, 
ptcxxutc  is  specified  on  an  unknown  body.  \o  accommodate 
this  boundarv  condition,  we  obtain  I  tom  Iteinoiilh's  eqna 
lion  and  the  peifecl  gas  law 


exp |  |  (  y  I  )■  y  I  ( v  <  !../> )  | 
l'l(-)  II  2|  A/,;  explKy  I )  y  |  ( s  *  | 

Since  the  entropy  v  is  constant  along  streamlines  and  is 
known  behind  the  tip  shock,  the  right-hand  side  ol  l  q.  (II) 
is  determined.  Therefore,  instead  of  l  q.  (9),  we  have  I  qs 
(4)  and  (II)  on  the  body,  lhe  flow  above  the  body  can  he 
calculated  exactly  as  before.  Hence  the  inverse  problem 
becomes  a  direct  problem  bv  the  method  presented  Ireie 

Numerical  Procedure 

Near  die  tip  of  the  body,  the  flow  is  taken  to  be  flow  over 
a  cone,  lor  the  rest  of  the  flow,  we  use  a  marching  scheme 
along  each  column  of  li.  A  good  approximation  is  first  ob 
tamed  at  the  points  on  that  column,  the  governing  equations 
are  then  iterated  until  an  error  tolerance  is  satisfied  before 
we  proceed  to  the  next  column 

It  is  worth  noting  that  we  are  tree  to  choose  anv  mesh  sia- 
of  d,  even  when  the  ("  and  (  characlenstics  are  nearlv 
parallel  Due  to  lq  (10).  we  are  also  not  loiced  to  lake  vetv 
small  mesh  si/es  in  o  either 

Approximalion  Solution 

I  gee  i  s  and  Savin1'1  noted  that  Its  pci  sonic  flow  ovci  axvm 
met  uc  bodies  can  be  approximated  as  locallv  two 


lleie  O  0  is  an  appioxiination  voiisistenl  with  the  approx¬ 
imation  R  0  I  liese  are  the  two  ordin.uy  diltcicntial  eqtia 
nous  which  give  us  the  iniital  guesses  in  the  itetation  scheme. 

Results  and  Discussions 

(  alciilaiious  ol  flow  ovci  bodies  ol  vanotix  shapes  for  a 
lange  ol  Mach  numbers  have  been  pci  formed.  The  average 
number  of  iterations  (geneiallv  between  2  and  4  tot  an  error 
loleinnce  ol  It)  ‘)  iisevi  pci  poim  decreases  as  the  total 
mi  mixer  ot  points  used  increases  In  view  ol  om  special  coor¬ 
dinate  system,  tew  points  aie  requited  to  describe  the  entire 
llowlteld  kind  points  ate  spaced  appropriately  according  to 
the  natural  sanations  of  the  flow  I  Ins  is  demonstrated  in 
lie  1  uiih  a  AO" « -thick  parabolic  body  ai  A/,,  2.  where  we 

compare  a  calculation  using  248  point'  in  the  whole  flow-field 
to  one  using  829  points. 

I  tguie  2  shows  the  piexxmc  dixti ibutioiix  on  two  parabolic 
bodies  I  lie  llow  neai  the  tip  is  conical,  hence  pressure  is 
constant  ihetc  lhe  presence  ol  a  piexxiirc  minimum  on  the 
body  should  tie  noted.  No  such  minimum  appears  in  two- 
dimensional  parabolic  wings.  Cousistem  xvuli  tlus  is  lhe  facl 
l liar  lhe  minimum  moves  further  down  along  the  body  as  A/„ 
oi  the  thickness  incieases.  because  thicker  bodies  arc  more 
two  dimensional  I  Ins  should  be  of  some  interest  since  the 
piesence  ol  a  pressing  minimum  in  some  cases  can  be  an  in¬ 
dication  of  flow  separation. 

I  igutc  f  shows  the  changes  of  R  (normalized  to  its  value 
ai  the  shock)  along  stteamhne  for  a  I <)“«  body  at  A/„  2. 

I  lie  streamlines  shown  lie  on  the  body  and  at  about  0. 1 . 0.5. 
and  I  body  length  away  I  torn  the  axis  As  one  can  see,  R 
changes  inoiiotonicallv  along  the  hodv  lhe  appioxiination 
R  constant  along  the  whole  both  is  pooi ,  but  R  does 
become  neailv  constant  at  less  than  a  body  length  away. 

lo  test  the  method  loi  solving  the  invoice  problem,  a 
direct  problem  was  solved,  and  the  resulting  pressure 
distribution  on  the  both  was  used  in  the  invoice  method.  I  Tie 
agreement  was  excellent  and  the  computational  time  is  com¬ 
parable  to  the  diievi  problem 

(  (inclusion 

N  stieamhnes  eli.ii ,n  ter isi ics  cooidmatc  system  for  ax 
isvnimelric  llow  lias  been  piecemed  In  om  coordinate 
xv stem  the  hodv  lies  along  one  coordinate,  and  the  shock  is  a 
straight  line  with  gml  points  tailing  exactly  on  it  Ranking 
Hiigonioi  conditions  aie  used  at  the  shock  C  onsequcntly , 
the  scheme  is  utheicnilv  accurate  Relatively  tew  points  are 
requited  lo  desiiihc  the  entire  I  low  liold  since  glut  points  are 
spaced  wccordtnr  lo  the  naluial  vaiiatioli  ol  the  tlow,  I  von 
in  legions  nt  neai  souloxccncc  ol  the  (  and  (  dial 

.uieiisiics.  we  aie  not  toued  lo  lake  veiv  sma'I  n.esli  si/es 
81.11  ting  liom  a  simple  iiiinal  guess,  ttie  scheme  vonveiges 
i  atiutlv  lo  lhe  exact  solution  Ihuauxe  so  lew  poinls  and 
iicialioiis  aie  needed,  [lie  method  is  sompiuuiioiiudv  veiv  el 
Inieni  (  )ui  hodv  til  i  ooidiu.Uc  svsiein  also  allows  us  lo 
solve  the  invoice  design  piohlcm  with  ease,  and  due  to  out 
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coordinate  system  the  inverse  problem  becomes  a  direct  pro¬ 
blem.  Our  results  for  thin  bodies  at  low  Mach  numbers  show 
there  is  a  pressure  minimum  on  the  bodies,  which  can  imply 
flow  separation  in  some  cases. 
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Supersonic  inviscid  (low  over  nonaxisynimctric  bodies  is  considered.  A  new  version  of  the 
method  of  reference  planes  is  used  In  this  version,  a  near  characterislics-streamlincs  coor¬ 
dinate  system  and  a  highly  efficient  numerical  integration  scheme  is  developed  The  CEL  con¬ 
dition  is  rigorously  satisfied  on  the  (low  Several  sample  calculations  arc  presented.  <  ivk7 

Academic  Press.  Inc 


1.  Introduction 

The  method  of  characteristics  for  three-dimensional  flows  has  been  developed  in 
a  number  of  ways.  Surveys  of  this  method  have  been  given  [1,2,3],  and  the 
leading  approaches  have  been  compared  [4],  The  main  advantages  of  such 
methods  lie  in  their  intrinsic  use  of  characteristics  as  well  as  their  accurate 
calculation  of  shockwaves.  Generally  these  methods,  which  require  consideration  of 
characteristic  conoids  and  bicharacteristics,  arc  regarded  as  complex  and  com¬ 
putationally  inefficient  compared  to  the  more  popular  finite  difference  shock  captur¬ 
ing  and  shock  fitting  methods,  c.g„  [5,  6,  7], 

Another  class  of  schemes  allied  to  the  characteristics  method  but  much  simpler  to 
apply  is  gencrically  referred  to  as  reference  plane  methods  [8-12,1].  Another 
designation  is  method  of  near  characteristics,  a  terminology  which  reflects  the  idea 
that  characteristics  arc  employed  in  an  approximate  fashion.  In  this  paper  we  apply 


*  Work  supported  in  part  by  grants  from  the  National  Science  Foundation  (C'UF  XI  04021  |  and  the 
ITS.  Department  of  Energy  (1)1:  ICOS  X5FR25000I 

’  Work  supported  bv  the  Air  force  Office  of  Seicntifn  Research  lAEOSR  5  2X3201 

378 

002I-‘W9I  ,  X7  $V00 

(  opytijihl  i  IVX7  b>  Aca<lcmi«.  Press.  Iik 
All  rights  of  reproduction  in  ans  form  rcsruftl 


,-Ti  -  A  f.m  f  *.  J.-X  A, 


...u/VT-vV-V 


,w  m  m  o ■  Pi  y  w  i  n  >u  u  ■  j  »j  pjt j  m  wm  v  «-  v  ■»-'.  -  -'* 


si  i’i  RsosK  iNvisnn  now  379 

a  variant  of  this  approach  to  the  problem  of  flow  past  nonaxisymmetric  bodies.  Our 
approach  is  most  closely  related  to  that  of  Sauer  [9]  and  Rakich  (12]. 

For  the  case  treated  here,  flow  past  a  body  is  divided  into  a  set  of  azimuthal 
planes.  In  each  plane  a  highly  successful  two-dimensional  characteristics  method 
[13,14]  is  applied.  The  "cross-talk"  between  such  planes  created  by  azimuthal 
derivatives  and  velocities  then  serve  as  forcing  terms  in  the  equations.  Unlike  earlier 
treatments  that  we  are  familiar  with  we  are  able  to  rigorously  consider  domains  of 
dependence  follow  the  Courant  Friedrichs  Lewv  ((’FT.)  condition  The  result  is  a 
method  which  is  extremely  fast  without  loss  of  efficiency  or  accuracy. 


2.  Formulation 

Since  the  form  of  (he  governing  equations  is  not  standard,  we  now  outline  their 
development.  Flow  in  cylindrical  coordinates  (  v.  /-.  tfi)  is  governed  by  the  following 
equations: 
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In  addition  to  the  continuity  (I  )  and  momentum  equations  (2),  (3),  |4).  we  have 
Bernoulli's  relation 


ir  + 1"  +  w  a*  M;,  1 

- + - =  ■ — -  + - .  ( 5 ) 

i  —  1  ->  -  _  1 

The  gas  is  specified  bv  the  state  equation 

—  <■  '  "v  =  constant,  (61 

l> 

where  the  entropy  S  satisfies 

u  VS  =  0.  1 7  ) 

between  shocks  while  (5)  is  also  valid  across  shocks  [  13  ].  We  normalize  u.  r. 
and  the  speed  of  sound  a  by  the  upstream  speed  of  sound  u„;  v,  r  by  the  body 
length;  p  by  its  upstream  value  /»„;  />  In  yp„:  and  .S'  is  replaced  by  (S'  S',, )  R.  where 
R  is  the  universal  gas  constant. 
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We  introduce 


W(.v,  r,  </>)-  tan 


the  flow  deflection  angle  in  the  projected  plane,  ^  =  constant.  Similarly 


ft(.x,r,</>)  =  sin  '  — — 


te) 


is  the  projected  flow  Mach  angle,  where 


,r  +  r‘ 


In  Appendix  A  it  is  shown 


sin  2/i  (  J ,  .S'  </.  HA 
</.  (»±  /*(/!))=  ±  — ^ ^ - A— f  ) 

G',  +  (tan  I)  tan  //  +  G':)  tl .  r 
tan  0  ±  tan  //  r 


C  *  :  <1 ,  =  ^-  +  tan(d  ±  //)  4- 

r.v  IT 


denotes  differentiation  in  the  C '  -directions,  i.c.,  in  the  near  characteristic  direc¬ 
tions.  ()±P(u\  are  the  corresponding  “Ricmann  invariants."  The  various  terms 
appearing  in  (11)  are  defined  in  Appendix  A.  In  what  follows  we  also  use 

”-l"',=7rTl„(l  +Lr  "*)/(' 

where  to  =  wiq.  The  second  form  for  a  in  (13)  follows  from  (5)  and  (6).  It  should  be 
noted  that  for  axisymmctric  flow,  we  have  w  =  0^  =  a#  =  0,  and  (II)  reduces  to  the 
appropriate  axisymmctric  equations. 

Next  we  define  new  coordinates  (a,  /i,  y)  such  that 

r,  =  v,  tan(d  +  ft ),  r/t  -  .v„  tan  (KZ  =  <j>.  (14) 


T  his  mapping  is  further  lixed  by  the  condition  at  the  body 


St  1*1  RSONK  INVISCII)  I  LOW 


.tin.)  the  condition  that  the  shock  have  unit  slope  at  each  constant  2  plane,  i.e.. 


at  the  shock  l ,'nder  this  transformation  the  (' 1  form  of  (111  becomes 


sin  2fi  j  .S',  U  ,  ,  -(!an(ltann+(/.|r. 


tan  0  +  tan  /r  r 


while  the  (  form 


sm  2//  ,  />.S  /)IC  \  O',  +  tan  0  tan  /i  -t  (P  Dr 

2  '  ;■  ;  I  *  tan  0  tan  //  r 


where  dilTerentiation  in  ('  direction  is  now  given  bv. 


r  2  tan  0  r,  r 

D  =--  - - — - -  — 

t  i  tan  II  +  tan  //  r/f  <  p 


As  shown  in  Appendix  B  the  4>  and  2  derivatives  are  related  by 

t  (r  tan  II  .v- )  L5(l  4  (tan(((  +  ft)  v  -  /•  )  v,(r  f/t)  r 
i'll i  (land  tanld  +  //))  v,  v/( 


If  we  combine  (17)  and  (IX).  then 


,  v  sin  2//  .S',  If, 

«,,=-■(!  tan  0  tan  /<  I  ~ — - —  — — --M/-(/t)l, 


.  ,  Sl”  I  11  ,<  V\ 


2  ("  22c!  „  _L  i„  ( u  t'-_l  ) 

■  sm  2/t  vi'  2  sin  u  I 


sm  2 fi  V 


_  sin  fi 


I  lie  entropy  equation  as  shown  in  Appendix  A  is  now 


S  —  -  .S' .  v,, , 
i  cos  0 


l 
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while  the  <j>  component  of  the  momentum  equation  is 


sin2  cos  0  u) 

- r r +  —  w„  +  -  ID* 


o)  sin2  it  cos  0  .  ,  v>  . 

■ - o  ,  -  ( I  +  a)  )  —  sin  0. 

A/,  /• 


With  the  use  of  (13)  a  can  be  eliminated  from  ( 22 )  to  give 


oj/(  =  to  tan  n  P ^ 


■  cos  (I 


(<u  sin  0  +  tail  /i  Pf) 


tan-  fi 


I  +  tan2  n 


r  cos  0  \  V  V  I  /  \  ;■  I  V  > 


■  cost) 


—  u>(o)j  +  or  sin  0 1. 


[231 


The  dependent  variables  v,  r,  0,  ft,  o>,  and  .8  are  determined  by  (14).  (17).  (IK)  or 
(20),  (21 ).  and  (23).  On  the  body  a ,r,0  are  specified  by  the  boundary  conditions 


v(  a  =  0, /u;)  =  //. 

/( x  =  0,  fi.  0=  /(/?.  2). 


(24) 

(25) 


<>[*  =  (),/;.:)  =  tan  1  ( ^  (/(.:> 
</> 


(26) 


In  addition  to  the  shock  conditions,  we  also  apply  (16)  at  the  shock. 


3.  Numerical  Pkoc  i  durfs 

The  following  scheme  is  an  extension  of  the  one  used  in  the  axisymmelric  case 
[13],  Each  azimuthal  section  (y*  =  constant:  A:  =  1,2,...)  has  the  (x. /()  grid  shown 
in  Fig.  I.  In  the  neighborhood  of  the  tip.  0  $/($/(,.  the  flow  is  taken  as  flow  past  a 
cone,  not  necessarily  circular.  For  /( >/(,  a  marching  scheme  described  next  is  used. 

Regard  the  flow'  as  determined  for  all  columns  up  to  /!  =  /(„  for  all  Ct  sections.  We 
first  indicate  how  the  flow  is  determined  at  the  body  of  the  /(,..  column,  denoted  by 
a |  in  Fig.  I.  When  integrating  in  the  C  direction  as  it  explicitly  appears  in  (18),  we 
trace  this  near  characteristic  back  to  the  /f„  ,  column,  the  point  /»,  in  Fig  1  Flow 
variables  at  hx  are  found  by  a  second  order  accurate  interpolation  scheme.  2 
derivatives  are  taken  as  the  averages  of  the  center  differences  at  </,  and  at  />,.  All 
other  equations  arc  integrated  by  appropriate  elementary  grid  points  differences 
since  they  only  involve  a  and  /(  derivatives.  The  values  of  the  flow  for  the  entire  row 
CA,  k  =  I....,  of  body  points  arc  now  iterated  until  a  convergence  criterion  is  met. 
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Above  the  body,  instead  of  the  boundary  conditions  (14)  (26)  we  have  ( 14 ),  (17). 
Again  values  for  an  entire  row  are  iterated  together.  We  proceed  in  this  manner 
until  the  point  a,„  of  Fig.  I  is  reached.  At  this  point  the  C  near  characteristic 
strikes  the  shock  before  hitting  the  //„  ,  column.  All  remaining  rows  in  /(„  are  now 
iterated  simultaneously. 

To  consider  stability  denote  the  true  Mach  angle  by  ji.  It  then  follows  that 

,  a 2  sin’ n 
sm  -  //  =  — - ^  . 

(/  +  If  I  +  <!)- 


and  hence 


ft  $  ft- 

At  each  constant  £  plane,  the  true  domain  of  dependence  for  the  flow  projected 
onto  that  plane  lies  inside  the  domain  of  dependence  determined  by  the  two  near 
characteristics.  As  for  the  £  derivatives,  let  b  in  Fig.  2  be  the  point  to  be  considered, 
and  let  4i  and  £:  be  the  angles  (in  ( v.  r,  <j>)  space)  between  r„,.  and  r,„  and  between 
r, and  rM.  We  want  £^ft  and  c:>  fi.  This  requirement  is  easily  satisfied  unless 
the  aspect  ratio  (the  ratio  of  the  largest  to  the  smallest  radius  at  the  cross  section)  is 
large,  in  which  case  a  smaller  step  si/e  of  / i  is  required 


t  i:,  2.  (((,  C)  coordinates  with  correspondinj!  angular  value  in  ( v.  r.  space 
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Sonic  mention  of  the  calculation  over  a  noncircular  cone  is  ordered.  Near  the  tip 
(0  ^  //  ^  /(,,  I  at  each  section  2  ~2k,  we  approximate  the  flow  as  being  the  flow  over 
a  circular  cone  with  half  angle  (),  =  tan  '(  /  t/f  =  /<,,,  C  =  st  1 1  We  then  compute  the 
rest  of  the  flow  (/(>/(,,  )  by  integrating  the  (nonaxisymmtric)  equations  using  the 
numerical  scheme  described  above.  Note  that  flow  over  a  noncircular  cone  is  con¬ 
stant  along  the  shock  and  the  body  at  each  section  C  =  Ca -  Hence  we  compute  the 
flow  along  the  /I  direction  until  this  constancy  condition  is  met  within  a  prescribed 
error  tolerance. 

Since  many  of  the  steps  in  our  procedure  are  iterative,  a  good  first  approximation 
can  significantly  accelerate  convergence.  To  motivate  our  choice  of  a  first 
approximation  observe  that  (17),  (18).  and  (21)  differ  from  axisymmelric  flow  in 
the  “second  order"  terms, 

,  r  cu) 

cr.  . 

ri/>  rtjt 

where  as  (19)  indicates,  f  <~(j>  =  OIck'c ).  Thus  if  a  body  can  be  regarded  as  locally 
axisymmelric  taking  the  flow  as  axisynimctric  should  be  an  excellent 
approximation.  In  any  case  locally  axisymmelric  flow  is  the  first  approximation 
adopted  by  us. 


4.  Ri  si  t  ts  and  Discussions 

For  purposes  of  exposition  we  have  performed  calculations  of  flow  over  bodies 
with  elliptical  cross  sections  and  azimuthal  parabolic  profiles.  Relative  few  points, 
planes,  and  iterations  are  needed  to  compute  the  entire  flow  field.  tf>  =  0,  n  and 
<p  =  ni2.  In  2  are  planes  of  symmetry  which  are  n  •(  assumed  in  the  calculations. 
and  arc  used  as  a  check  on  the  correctness  of  the  results.  Throughout  the  entire 
flow  field,  all  grid  points  with  such  symmetry  are  found  to  have  values  in  agreement 
within  the  same  order  of  magnitude  as  the  prescribed  error  tolerance. 

Figure  3  shows  the  body  and  shock  along  the  half  planes  <f>  =  0  and  <j>  =  n  for  flow 
at  A/,,  =  2  over  a  body  with  30%  thickness  at  <j>  =0  and  20%  thickness  at  <f>  =  n  2. 
To  carry  out  this  calculation  we  took  32  a/imuthal  sections  each  having  164  grid 
points,  20  on  the  body.  The  result  of  reducing  the  number  of  azimuthal  planes  to  16 
is  shown  by  v  s  and  +  "s.  For  calculations  w  ith  the  fine  mesh  size,  the  lines  for  </>  =  () 
and  <t>  =  n  are  indistinguishable  in  this  figure.  For  calculations  with  the  coarse  mesh 
size,  +  signifies  a  point  at  </>-()  and  x  signifies  a  point  at  </>  =  n.  Fach  pair  of  + 
and  x  appearing  together  in  the  figure  have  the  same  values  of  x  and  ft.  and  hence 
they  have  the  same  flow  values.  Figure  4  shows  the  body  and  the  shock  for  the 
same  flow  along  <j>  —  7t  4  and  </»  =  3 it  4.  Again,  for  calculations  with  the  fine  mesh 
size,  the  lines  for  the  two  half  planes  are  indistinguishable.  For  calculations  with  the 
coarse  mesh  size,  +  signifies  a  point  at  <j>  =  ~  4.  and  x  signifies  a  point  at  <f>  -  3n  4. 
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with  each  pair  of  +  and  x  having  t he  same  i  and  //  In  these  two  figures, 
agreement  is  excellent  with  respect  to  both  symmetry  and  the  different  grid  si/es. 

For  the  same  flow,  the  pressure  distribution  on  the  body  at  <t>  (I.  n  4.  n.  4  is 
shown  in  Fig.  5.  For  <f>  ~  0.  n  the  agreement  is  excellent  I  -  <sr  </>  =  n  4.  .Ft  4  the 
agreement  is  good  with  respect  to  symmetry,  but  there  is  a  small  discrepancy 
between  the  coarse  and  the  fine  grid  calculations  Since  the  grid  si/e  in  the  direc¬ 
tion  is  quite  large,  about  rt/O,  the  discrepancy  is  certainly  tolerable. 

The  cross  flow  on  the  body  is  given  in  Fig.  6  As  the  figure  shows,  the  cross  flow- 
vanishes  in  the  symmetry  plane.  At  <j>  -  n  4.  3 r.  4.  there  is  some  discrepancy  between 
the  coarse  and  fine  grid  calculations  due  to  the  largeness  of  the  C  grid  si/e.  Figure  7 
shows  that  cross  flow  at  the  shock  The  results  are  again  quite  good.  Figures  X  and 
9  show  the  same  body  at  two  degree  angle  of  attack  traveling  .it  A/,,  =  2 

As  an  indication  of  the  computational  speed,  we  mention  that  for  the  whose  flow 
field,  the  time  of  a  typical  calculation  is  roughly  It)  sec  on  an  IBM  30X1. 


Al'iMNDix  A  Fquaiions  in  Nonaxisymmi  ikk  Flow 

We  wish  to  express  the  governing  equations  tit  |7|  in  near  characteristics  form 
Since  p  =  p( p.  .V).  using  (6 1  and  (7).  ( I  )  becomes 
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Denote 


q  =  (i/.t)  and  V\  =  (4-.4r) 

\r\  i  r ) 

I  hen  (2)  and  (3  I  become 


q  ■  V \q  +  -  V.  />  =  f. 


II  ('  /  I!  ‘  \ 

In  a  constant  half-plane  with  0  and  /i  defined  by  (X)  and  (9).  let  t((l)  be  the 
tangent,  and  n(//)  be  the  normal  on  the  projected  streamline  Denote  v  and  n  as  the 
arc  length  and  normal  coordinate  I  hen  <  A  I  )  becomes 

<P  ( >' ‘I  <'»  \  {  "  <>  .  ;i”  ,  IP \ 


p  “i  i">n  .w  ir 1 "  \ 

—  +-  :P  -  +  </--)  -  -  -  t:  +  —  +  ^7 

i  v  V  <  v  (  n )  \  /•  (  q>  r  r  <  </>  > 


and  (A2)  becomes 


('<;  I  ('/>  -  it  rii  n  . 

(/  f  -  —  - -7—  4  —  sin  ll. 

<m  /i  (  \  r  i  </>  r 

,  i'll  1  i'p  -  u</  ('ll  u 

= - —  f  —  cos  0. 

is  1 11  11  /  cif)  1 

Denote  <r  =  In/)  and  id  =  it  </  Rearranging  (A3).  (A4|.  and  (A5). 


(<» ,  ' 

i'll  i  ,  (A/',  -  1  )' 

'  I  1(7  1 

v  <  v  -  <  a/  ;  i)1’ 

'  r.v  "  (A/;  -  1 1 

\ 1  •  rn  ) 

ll1 

where  /•',  -- (o'  1 )  cos  II  {u>,r)tl,. 


I  id  sin  0  , 

l-  •  --  -  id  t  —  a  t  + - (I  for) 

r  '  r  r 


(A6)  can  then  be  written  as 


sin  2/i  /',  cos  n  +  /  .  sin  11 

il  ,  II  ±  — - - ll  .  ci  - - : - - - il , 

2y  sin(fl±pl 


where  </ ,  is  given  in  (12 ). 
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Now  we  give  a  similar  treatment  lor  the  <t>  momentum  equation  (4).  This  can  be 
written  as 

v/>  I  i  n  u  i  n  m 

(q  ■  V, )  n-  + - —  ■--- - — - . 

p  )rp  i  <p  r  rip  r 

In  natural  coordinates,  this  becomes 

i  (/(■)  if  in  qu)  rift  it  ,  <n  sin  <p 

q  — - 1-  —  — - — ~  q - 

(  s  -.  /•  r  i  ip  r 

Using  (A4)  and  (AI3),  this  then  becomes  (22). 


Aimmndix  B.  Coordinati  Trassi-ormatiom 
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< M  y.  /I. 

2), 
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the  chain  rule  we  obtain 
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i  /  ’• 

y. 

\ 

I*. 

p. 

/ 

\  it> ,  </),,  d  1 

1  \c. 

s  r 

<!<  1 

/ 

Since  </S  =  s",  and  <fi,  --=  <j>v  =  0,  <!>  =  I. 


where 


Then 
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Steady,  invisod,  irrolational  firm  of  a  perfect  pas  in  two  dimensions  is  considered  in  the 
tangent  gas  approximation  A  fast  and  accurate  method  of  solution  is  proposed  and  solved 
numerically.  Comparison  of  tangent  gas  and  exact  flows  are  presented.  Tangent  gas  solutions 
when  used  as  the  first  step  in  the  iterative  solution  of  the  exact  flowfield  are  shown  to  give 
substantial  reduction  in  computational  lime.  1  ivxs  AuUcmic  Prr>».  In, 

1.  Introduction 

The  computation  of  steady  flow  past  an  airfoil  is  crucial  to  the  determination  of 
aerodynamic  characteristics  such  as  lift,  drag  and  moment  coefficients.  In  many 
instances  potential  theory  suffices.  Neglecting  viscosity  it  is  exact  for  shockless  flow 
and  is  a  satisfactory  approximation  for  transonic  flow  with  weak  shocks.  For  two 
dimensions  the  calculations  arc  usually  carried  out  in  a  conformally  mapped  plane, 
an  approach  used  by  Sells  [1],  Garabedian  and  Korn  [2],  and  Jameson  [3], 
Similar  techniques  have  been  used  for  multi-element  airfoils  [4,5]  and  nacelies 
[6]  Three  dimensional  potential  theory  has  been  treated  by  Caughcy  [7]. 

Since  the  equations  arc  nonlinear,  the  potential  equation  is  usually  solved 
iteratively.  In  some  instances  the  potential  equation  does  not  admit  unique 
solutions  [8  10]  and  in  addition  becomes  a  poor  approximation  for  increasingly 
strong  shock  strengths.  As  a  result  more  recent  investigations  treat  the  full  Euler 
equations.  Finite  difference  and  finite  volume  methods  have  been  successfully 
implemented  by  Jameson  [11]  and  l.crat  and  Sides'  [12],  Because  of  slow  rates  of 
convergence  considerable  effort  has  been  directed  towards  accelerating  these 
methods  [13].  Convergence  rates  depend  on  factors  such  as  the  grid,  initial  guess, 
time  stepping  scheme  and  method  of  solution. 

*  Prevent  address:  (intrant  Institute  <>t  Mathematical  Sciences.  251  Mercer  Street.  New  York. 
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In  this  paper  we  present  a  set  of  flow  dependent  grid  systems  and  initial  flowfield 
guesses  which  substantially  improve  convergence  rates  when  applied  to  the  Euler 
equations  for  flows  past  an  airfoil.  These  are  based  on  solution  of  the  tangent  gas 
equations  introduced  by  Chaplygin  [14]  and  further  developed  by  von  Karman 
and  Tsien  [  15,  16]. 

Woods  [17],  who  extensively  studied  these  equations,  proposed  certain  iterative 
methods  for  solving  both  the  analysis  and  design  problems  for  flows  past  an  airfoil. 
The  methods  developed  in  this  paper  are  substantially  different  and  offer  a  method 
for  a  fast  and  accurate  solution  to  a  problem.  (We  have  also  addressed  the  inverse 
problem  and  presented  an  exact  method  for  its  solution  [18].) 

As  will  be  seen  the  tangent  gas  solution  lies  close  to  the  Euler  solution  even  for 
high  subcritical  flows.  This  is  used  as  a  basis  for  iterative  solution  of  Euler  equation 
for  flows  past  an  airfoil  by  means  of  FL052S  (written  by  A.  Jameson,  E.  Turkel 
and  M.  Salas).  The  grid  used  is  the  natural  one  generated  by  the  tangent  gas 
equations  and  the  starting  guess  is  the  tangent  gas  solution.  As  will  be  seen  this 
results  in  substantial  computational  reduction  even  for  supercritical  flows. 


2.  Basic  Equations 

Consider  steady,  inviscid.  irrotational  How  of  a  perfect  gas  in  two  dimensions, 
then  in  the  usual  notation 

V  ( i></ )  •-  0.  V  x  q  =  0.  pif>‘  =  I .  ( 1  ) 

The  variables  are  normalized  by  their  free  stream  values  and  linear  dimensions  by 
an  appropriate  lengthscalc. 

The  stream  function  ip  and  potential  4‘  ;,rc  introduced  in  the  usual  way 

/K/  =  rVx(i/ik),  ^=V<*,  (2) 


where  k  denotes  a  vector  perpendicular  to  the  plane  of  motion.  The  constant  c  has 
been  introduced  for  later  purposes. 

If  .v  and  n  are  local  distances  along  streamlines  and  potential  lines,  respectively, 
(2)  can  be  written  as 


ds  -t-  i  dn  =  -  ( dip  +  i  —  dip 

‘i  \  n 


If  equations  can  be  derived  that  map  the  space  of  <j>,  ip  on  to  the  space  of  the 
velocity  magnitude  and  direction  (</.  0),  then  one  can  take  advantage  of  the  fact 
that  the  tangent  of  the  flow  direction,  tan  0.  is  the  same  as  the  slope  of  the  airfoil 
surface  where  i //  -  0.  Then  if  q  vs.  0  can  be  found  corresponding  to  i //  =0  on  <f>,  i p 
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plane,  the  state  of  flow  on  the  airfoil  surface  will  be  known.  Toward  this  end,  we 
write  Eq.  (3)  alternatively  as 

HI  /  \ 

dz  —  dx  +  /  dv  =  —  (  d<t>  +  i  -  dip  ),  ( 4  ) 

<i  \  P  1 

where  v,  y  are  cartesian  coordinates  and  0  flow  direction  angle.  With  q  and  I)  as 
independent  variables,  it  is  easy  to  derive  from  (4) 

,  q  ,  ,  !  -  M1  , 

<t>o  =  ~  <’V«.  - —  c4>„.  (5) 

p  pq 

If  dependent  are!  independent  variables  are  interchanged  and  the  Prandtl  Meyer 
function 


is  introduced  in  place  of  q,  then 


>V=0,  0^  ±  A'(v)  v*  =  0. 


The  +  sign  refers  to  subsonic  and  supersonic  conditions,  respectively  and 


All)*/! 


/>(</(  A/))' 


I  epical  physical  ;(  -  v  -t  iy)  and  potential  ir(  =  <f>  +  itp  )  planes  are  shown  in  Fig.  1 . 
I  he  airfoil  maps  into  a  slit  in  the  u-plane.  The  gap  Hli'  in  the  potential  plane 
corresponds  to  / '.  where  circulation  about  the  airfoil  is  -  /'. 

I  he  system  ( 7 1  should  be  solved  subject  to  the  density  speed  relation  obtained 
Irom  1 1  l  and  Bernoulli's  relation 

q  1  r  <//> 

—  +  —  —  I  —  =  constant.  (10) 

: •  p 


3.  I  AM.IM  C.AS  AI'I'RoXIMaTION 


I  v|uations  1 7 1  are  nonlmeai  and  are  therefore  diflicult  to  .  've.  A  good 
appt ■  'Viinal ion  to  those  equations  under  certain  conditions  can  be  obtained  bv 
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PlO  1  Airfoil  in  physical  r-plane  and  potential  it-planc. 


introducing  the  so-called  “tangent  gas  approximation”  [17],  in  which  the  iscntropic 
relation  between  p  and  p  given  in  Eq.  (1 )  is  replaced  by  a  tangent  to  the  curve  of  /> 
vs.  1/p.  This  approximation  is  then  given  by 


(ID 


From  ( 10)  we  obtain 


/’  =  /<//(  , 


(12) 


With  the  constant  c  in  (8)  taken  as 


<  =  1/0.  , 


we  obtain  from  (8) 


A  ( v )  =  I. 


03) 


(14  i 


Then  for  subsonic  flow  (7)  becomes  the  Cauchy  Picmann  equations 

»,  +  *•,  =  (>.  115) 

Equations  (15)  arc  exact  for  the  tangent  gas  and  also  for  incompressible  flow 
(M  =  0).  In  addition,  it  will  be  seen  that  it  can  be  a  very  good  approximation  to  the 
original  equations.  In  the  above  formulation  the  langency  point  has  been  taken  to 
the  freestream 
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With  this  selection  of  tangency  point  (he  following  relations  hold  for  the  tangent 
gas  [17] 

q  ~  sinh  v *  cosechfv*  v),  /<  -  tanh(v*  v|,  t  =  - - - - — .  (17) 

I  /(,  colh  v 

where  the  contant  v*  is  given  by 

'■''"(rr;) 

From  (6)  it  is  seen  that  v  ,  =  0  and  at  stagnation  points  (denoted  by  zero  subscript  I 


4  S<lf  IT  ION  PR(K  iim  Rt 

It  follows  from  (15l  that 

r  =  v  +  ill  ( 20 ) 

is  an  analytic  function  of  u  It  will  be  useful  to  map  the  vc(  =  <t>  +  up  I  plane  onto  the 
plane  of  a  new  variable  a  =  |<r|  en  such  that  the  body  in  the  w-plane  which  is  a  slit 
(a  part  of  the  line  tp  =  0)  maps  onto  the  unit  circle  n  =  c";  0  $  or  ^  2n  and  the  rest  of 
the  w-plane  maps  onto  the  exterior  of  the  unit  circle  This  is  accomplished  by 

w  =  a{ae  +  a  V’">  +  i2a  sin  a„  ln(<re  '*“)  (21) 

which  allows  for  angle  of  attack  and  circulation  about  an  airfoil  surface,  to  be 
related  to  |<r|  =  1.  Circulation  - is  related  to  the  constant  a  by 

r  =  4n«  sin  otn.  (22) 

Here  constants  a  and  a„  arc  as  yet  unknowns. 

From  (21 )  one  obtains 

^  =  I  -o  1  )(<■  '*•-<1  ')•  (23) 

tin 

On  the  body  a  =  c”;  0  s;  or  <  2 n,  <p  and  ip  arc  given  by 

^(or)  =  2r/[ cos(a  -  ot„ )  -  (a  -  or,,)  sin  or,,],  i p(i)  =  0  (24) 


or,  in  (23 )  is  given  by 


,»4  Vl  a* A  .*1  *«i 


,«l  >t|  .««  , 


,U«>  ■*!  J<  ij 
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Thus  the  rear  and  front  stagnation  points  map  into  a  =  i  and  a  =  e"\  respectively. 

Since  r  is  an  analytic  function  of  <r,  a  convenient  representation  of  r(a)  is  given 
by  (see  also  Ref.  [  19] ) 


cxp(r(<r))  =  ( I  —  fl  ')  "‘-o'  ')  '  exp  V  c„a  "l. 


where  <5  =  0,/n,  0,  the  trailing  edge  angle.  The  complex  constants  c„  arc  represented 

by, 

<'„  =  A„  + ili,,.  (27) 

Note  that  (26)  contains  the  Kutta  condition.  Two  Schwarz  Christoffel  factors 
appear  in  (26)  because  of  the  discontinuity  in  0  at  the  two  stagnation  points. 

From  (26)  the  relationship  between  upstream  flow  direction  0  ,  and  a0  is  given 
by 

0,  =B(l  +  7r  +  2x().  (28) 

The  free  stream  condition  is  given  by 


/«„  =  o. 


On  the  unit  circle,  (26)  reduces  to 


xp( T(t>'3' ) )  =  G(z)c"l,yl  exp  (  £ 


C/(a)=  2  sin-  |2(sin  a„  +  sin(x  -  a„))| 


'  =  ^(1  -<5)(rt-a)  +  ^a  +  ^  -  nV 


( a  —  a , )  +  at,, 


U(a  —  a,)  in  (32)  is  the  unit  step  function.  The  tangent  angle  0H  of  the  body  is 
related  to  0  by 

0(x)  =  0H(x)  —  n  -  ni/(x  -  2 , ).  (33) 

Separation  of  (30)  into  real  and  imaginary  parts  leads  to 

v(»)=  Y  M„eos  ii2  +  B„  sin  nx)  (34) 

n  0 

f!( 2 )  =  Y  ( Bn  cos  —  -4 „  sin  nx )  t  n  +  ( 35 ) 


r- 


pt  j 
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where 


and 


v(ot )  =  -  v(a)  -  In  O'(a), 


ff(a)-0„(cc)  -  -(!-<$)(/[-*)  -  (  a  +  ^ 


(36) 


(37) 


The  closure  condition  of  the  airfoil  is  related  to  the  leading  terms  of  the  series  by 
(Appendix  A), 


A  i  =  (I  -  S)  -  ( I  -  /( ,  )  2  sin'  a,, 
=(1  ~P,  )  sin  2a,,. 


(38) 

(39) 


5.  Analysis  (Diric  t)  Pkoiu  i  m 


Here  the  flow  past  an  airfoil  is  sought.  An  iterative  method  of  solution  similar  to 
the  one  for  incompressible  flow  (15)  is  found  to  converge  with  good  accuracy.  The 
method  of  solution  goes  as  follows. 

An  initial  estimate  of  arclcngth  as  a  function  of  circle  angle,  ,v(a),  (e.g„  of  a  flat 
plate  in  incompressible  flow)  is  made.  From  the  given  contour  0H(s).  ()fl(a)  is 
estimated  and  (7(a)  is  calculated  from  (37).  a0  is  obtained  from  (28).  After  the 
closure  conditions  (38)  and  (39)  are  imposed,  a  new  form  of  (7(a)  is  generated  and 
then  its  conjugate  v(a)  is  obtained  from  (34).  v(a)  is  then  obtained  from  (36)  and 
speed  q( a)  is  obtained  from  (17).  The  updated  value  of.v(a)  is  now  obtained  from 
q(a)  using  the  relation 


i  m 

q  da 


da 


=  2„f 


I  sin  a„  +  sin(a  —  a„)| 


da. 


(40) 


where  the  constant  a  is  now  given  by 

,v(2ji)  =  I.  (41) 

The  above  procedure  is  repealed  until  convergence  is  obtained.  The  criterion  for 
convergence  was  taken  to  be  that  maximum  difference  in  arc-length  between  suc¬ 
cessive  iterations  be  0(10  (’).  Typically  the  number  of  iterations  required  was  no 
more  than  eight  and  the  computation  time  was  roughly  one  second  on  an 
IBM  3081  with  128  points  taken  on  the  unit  circle.  The  actual  numerical  calculation 
is  facilitated  through  the  u->e  of  the  fas;  fourier  transform  (F'FT)  and  the  fact  that 
(34)  and  (35)  are  conjugate  fourier  series.  The  fourier  constants,  c„,  are  also 
obtained  easily  during  I  F!  which  are  used  for  generating  grids. 
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6.  Grid  Generation 

The  physical  plane  is  related  to  the  circle  plane  through  [17] 

‘'I  =  2C{(,+/'- 1-' £*-" -/'»<•  •%*]■  <«) 

Here  an  ovcrbar  denotes  complex  conjugate.  Note  that  for  incompressible  flow  r  is 
an  analytic  function  of  er,  as  it  should  be. 

From  (21)  and  (26)  it  is  easily  seen  that 


Pic;.  2.  Comparison  of  tangent  gas  solution  and  Pulcr  solution  over  NACA00I2  Airfoil  at 
Mach  =  0.6  and  angle  of  attack  0.0.  — ,  tangent  gas  solution:  +  +  +  ,  Ruler  solution 
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=  -ae  +  ,a°(  1 


(e  o  1 ):  exp  (  -  £  W 


Equations  (42),  (43)  and  (44)  are  used  to  map  the  circle  plane  into  physical  plane 
and  the  flowfield  variables  are  obtained  from  (26),  (17),  and  (18). 

Observe  that  the  grid  generated  is  flow  dependent.  Since  the  mapping  from  a 
plane  to  z-plane  is  not  conformal  except  when  M  —  0,  the  grid  generated  in  physical 
plane  is  not  in  general  orthogonal.  The  grid  produced  by  this  method  appears  to  be 
more  natural  than  the  incompressible  conformal  grid. 


o  0  2  0  4  0  6  0  8 


Kig.  3  Comparison  of  tangent  gas  solution  and  F.ulcr  solution  over  NACA  0012  Airfoil  at  Mach  0  7 
and  angle  of  attack  =0.0. - ,  tangent  gas  solution;  +  +  +  ,  Ruler  solution 
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7.  Rhsults 

Figures  2-5  compare  the  tangent  gas  solution  with  the  converged  Euler  solution 
(as  calculated  by  FL052S).  The  tangent  gas  solution  is  seen  to  be  remarkably 
accurate  even  at  the  near  critical  case  depicted  in  Fig.  3  and  the  slightly  critical  case 
shown  in  Fig.  4.  Even  when  a  clear  shock  is  present  as  in  Fig.  5,  the  tangent  gas 
solution  only  fails  in  a  relatively  small  neighborhood  of  the  shock. 

Figures  6  and  7  indicate  for  two  typical  cases  the  number  of  iterative  cycles  to 
achieve  a  convergence  criterion.  The  criterion  used  is  the  enthalpy  error  introduced 
by  Jameson  [20],  In  each  figure  we  indicate  the  number  iterations  required  to 
reach  the  indicated  criterion.  The  first  column  of  each  figure  refers  to  use  of  the 
tangent  gas  grid  and  the  tangent  gas  solution  as  a  starting  flow.  The  second  column 


0  OC  04  06  Ob 
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t-Ki  4  Comparison  of  tangent  gas  solution  and  I  nter  solution  over  NACAOOi:  Airfoil  a: 
Mach  -  0.50  and  angle  of  attack  5.0  degrees  .  tangent  gas  solution,  i  a  *  .  toiler  solution 
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gives  the  analogous  values  using  the  conventional  grid,  viz.,  that  generated  by  con¬ 
formal  mapping  and  a  uniform  flowfield  as  the  starting  guess.  (Little  change  in  con¬ 
vergence  was  observed  if  incompressible  flow  was  taken  as  the  initial  guess.)  As  is 
seen  the  reduction  in  cycles  is  substantial.  In  this  same  vein  if  the  convergence 
criterion  is  reduced  by  a  factor  of  10  the  comparison  becomes  more  dramatic  the 
tangent  gas  approach  leads  to  a  10-fold  reduction  in  cycles  over  the  usual  approach. 

In  order  to  distinguish  whether  the  grid  or  the  tangent  gas  approximation  was 
more  significant  in  speeding  convergence,  we  also  ran  the  programs  using  the 
tangent  gas  grid  with  a  uniform  first  guess.  Although  some  improvement  resulted, 
the  clear  implication  from  this  was  that  the  tangent  gas  solution  as  a  first  guess  was 
the  most  important  factor. 
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t  ic;.  5.  Comparison  of  tangent  gas  solution  and  fiulcr  solution  over  NACAOOIJ  Airfoil  at 

Maeh  -  0  758  and  angle  of  attack  —0.14  degrees  - ,  tangent  gas  solution;  (+,>')  fiulcr  solution.  +  . 

upper  surface;  ( ),  lower  surface. 


Fig  6  Euler  solution  (FL052S)  for  near  critical  now  past  an  NACA  0012  Airfoil  at  Mach  0  50  and 
angle  of  attack  -  5  0  degrees.  (  +  ,  O):  grid.  64«32;  grid  type,  tangent;  initial  guess,  tangent;  number  of 

cyclccs,  344  ( - );  grid,  64»32;  grid  type,  conformal;  initial  guess,  uniform;  number  of  cycles  913. 

Average  error  in  enthalpy,  0. 1385E-03.  +,  upper  surface;  O,  lower  surface. 


0  02  04  x  06  06  1 


Flo.  7.  Euler  solution  (FLOS2S)  for  supercritical  How  past  an  NACA  0012  Airfoil  at  Mach  0 .758  and 
angle  of  attack  =0  14  degrees.  (  +  ,  O):  grid,  64«32,  grid  type,  tangent;  initial  guess,  tangent;  number  of 

cycles,  381  ( - ):  grid,  64»32;  grid  type,  conformal;  initial  guess,  uniform;  number  of  cycles.  715 

Average  error  in  enthalpy,  0  2454E-03  +  ,  upper  surface;  C ',  lower  surface. 
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Appfndix  A:  Closurf.  Conditions 

If  C  is  a  closed  contour  around  an  airfoil,  the  the  closure  condition  is 

|  d:  =  0.  (Al) 

Hence  from  (42)  we  obtain  (see  also  Ref.  [17]) 

(I  +  0»)  |  da  =  (1  -fiy  )  |  e  (A2) 

From  (43)  and  (44)  it  follows 

e'  ~  =  -ae,m'  * [  1  +  —  +  0(o  2 ) 
da  l  ° 

c  ’  ~  =  W'*1  H  -e  2,’'  +  ^  +  O(0  2) 

da  L  a 

where 

K,  =  c,  +  <5  -  1,  K2  =  ( 1  +  S  +  r, )  e  ~  2‘*'  +  2e ~ (A5) 
Use  of  residue  theorem,  (A3)  and  (A4)  reduces  (A2)  to 

(I +/M**’*,  =0-0,)*  O0*2.  (A6) 

Equating  real  and  imaginary  parts  wc  obtain 

M,  +  <5  —  I ) cos  a0  =  B,  sin  a0,  (/l,  +  £  +  1  -  2(i,  )  sin  a0  =  fi,  cosa0.  (A7) 
From  (A7)  we  obtain 

A,  =(!  -<5(  —  (1  -//,  )  2  sin2  an,  5,  =  ( 1  -0,  )sin  2a0.  (A8) 

For  the  incompressible  case  (/(,  =  1 )  this  reduces  to 

A  ,  =  (!— <5),  fi,=0.  ( A9) 


(A3) 

(A4) 
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The  inverse  problem  in  the  tangent  gas  approximation  is  considered  An  exact  method  for 
designing  airfoils  is  presented.  Constraints  on  the  speed  distribution  are  easily  implemented.  A 
simple  numerical  algorithm  which  is  fast  and  accurate  is  presented  Comparison  of  designed 
airtoils  using  the  tangent  gas  method  with  exact  l-uler  results  is  found  to  be  excellent  for  sub- 
critical  flows  <  IVXh  Academic  Press.  Inc 


I.  Introduction 

As  is  well  known  (see  [I])  certain  types  of  pressure  distributions  achieve 
aerodynamically  desirable  features  such  as  delay  of  transition  and  boundary  layer 
control.  The  determination  of  an  unknown  airfoil  from  a  specified  pressure  dis¬ 
tribution  is  known  as  the  inverse  problem. 

Numerous  methods  for  the  two  dimensional  incompressible  case  exist  [2  5]. 
Compressible  inverse  methods  arc  for  the  most  part  based  on  some  kind  of  iterative 
procedure,  relying  on  either  a  Dirichlct  or  Neumann-type  boundary  condition.  In 
the  Dirichlct  formulation  [6-11]  a  sequence  of  boundary  value  problems  for  the 
velocity  potential,  with  wing  geometry  updated  at  each  step,  is  solved.  The  updated 
condition  arises  from  the  normal  velocity  resulting  at  each  unconvcrged  step.  For 
the  Neumann  formulation  [12  15]  a  sequence  of  analysis  problems  arc  solved  over 
a  corresponding  scries  of  geometries,  liach  geometry  is  provided  by  some  rational 
method  depending  on  the  error  being  driven  to  zero.  A  complete  survey  of  such 
methods  has  been  given  by  Slooff[!6]. 

In  this  paper  we  present  an  exact  melhod  for  two-dimensional  subsonic  flow 
within  the  limitations  of  the  tangent  gas  approximation  (17  19].  Woods  [20] 
extensively  studied  these  equations  and  proposed  certain  iterative  methods  for  solv¬ 
ing  both  the  analysis  and  inverse  problems  We  presented  a  substantially  different 

•  Present  address:  Courant  Inslilnlc  of  Malliemalic.il  Sciences.  2M  Mercer  street.  New  York, 
NY.  10012 
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method  for  the  anal' sis  problem  [21].  The  inverse  method  developed  here  is  non¬ 
iterative  and  exact. 

As  is  shown  in  [21  ]  the  tangent  gas  solution  lies  very  close  to  the  Kuler  solution 
even  for  high  subcritical  flows.  Therefore  the  design  of  an  airfoil  in  this  regime  by 
our  method  should  be  an  almost  correct  airfoil. 

In  this  paper,  we  have  been  able  to  show  that  the  direct  Eiulcr  solution  over  the 
designed  airfoil  is  very  close  to  the  input  speed  distribution.  Moreover,  the  con¬ 
straints  necessitated  by  upstream  condition  and  closure  requirements  are  very  easily 
incorporated. 


2.  If  ASK '  lot  A  I  l<  INS 

Consider  steady  two-dimensional  flow,  then  in  the  usual  notation 

V  |  /n|  i  0.  V  s  q  o,  p  p  I  ( I  | 

The  variables  are  normalized  by  their  free  stream  values  and  linear  dimensions  by 
an  appropriate  length  scale. 

The  stream  function  ip  and  potential  <!>  are  introduced  in  the  usual  way 

P q  <  V  <  1  i//k  |.  q  W.  1 2 1 

where  k  denotes  a  sector  perpendicular  to  the  plane  of  motion.  The  constant  <  has 
been  introduced  for  later  purposes. 

If  v  and  n  are  local  distances  along  streamlines  and  potential  lines,  respectively. 
(2  l  can  be  written  as 


</,s  +■  t  tin  -  -  ( </</>  +  i  -  dip  ) . 

</  V  1 


Alternately,  we  can  write 

d:  —  </.v  +  id\  -■  —  (  dip  +  /  -  dip  ) .  (4  I 

</  V  /’  / 

w  here  \  and  i  are  cartesian  coordinates  and  (>  the  flow  angle  If  </  and  0  are  taken 
as  independent  variables,  then  it  is  easy  to  derive  from  (4)  that 

,/  I  M:  , 

Ip, I  -  <  |//  <p,  - -  cip„.  (-■'  I 

r  /'</ 

If  dependent  and  independent  variables  are  interchanged  and  the  Prandtl  Meyer 
function 

pi  ,  ,  ,  dll 

v  |  (|l  A/"| )  •  -  lb) 
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I  ip  I  Airfoil  in  physical  -  plane  and  poicnlial  h -plane 


is  introduced  in  place  of  </.  then 


~ — i  'V  “  *)•  "..-f :*'<»•)  >>  =  ().  (7  | 

A  ( v ) 

The  t  sign  refers  to  subsonic  and  supersonic  conditions,  respectively,  and 


A't  >•  I  /> 


>Uf{  Mil' 


w  here 


/>••'  il  At 


(SI 


Typical  physical  r(  -  v  +  n  )  and  potential  u  (  < j>  +  np )  planes  are  shown  in  Fig.  1. 

The  airfoil  maps  into  a  slit  in  the  u  plane.  The  gap  HB'  in  the  potential  plane 
corresponds  to  /,  where  circulation  about  the  airfoil  is  /' 

The  system  (7)  is  augmented  by  the  density  speed  relation  obtained  from  ( I  )  and 
Bernoulli's  relation 


</  I  ;  <ip 

-I  —  -  —  constant. 

-  /  At',  r 


3  Maim’IM.s 


(  10 1 


The  u -plane  is  mapped  onto  the  exterior  ol  the  unit  circle  in  the  a  plane  by 

it  (itoc  ""  t  <7  1  i2ii  sin  jr„  ln(r7C  "‘  I.  (Ill 
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/'  is  related  to  £7(.v  =  1 )  and  Q(s  =  s, )  by 

/'=  2Q{s  =  sf)  -  Q{s  =  I )  =  4 na  sin  a0.  (21 ) 

Here  s,  denotes  the  distance  of  the  front  stagnation  point  from  the  upperside  of  the 
trailing  edge.  Q(s  -  I ),  (?(.r  =  .rf),  and  hence  /  are  known  from  the  given  surface 
speed  distribution  q,(. r). 

From  (20)  and  (21 ) 


Q{s=  I )  2 

- - —  =  -  (a„  +  cot  *„). 

/  n 


(22) 


i 


After  (22)  is  solved  for  a0  we  obtain  the  constant  a  from  (21 ).  Next  (5(x)  is  com¬ 
puted  from  (19)  and  v(a)  is  obtained  by  inverting  ( 18) 

*(ot>  =  0  ‘(0(a))  (23) 

and  </3t)  =  </J.v(a))  is  obtained  from  (17). 

Thus  far  our  deliberations  are  exact.  Ideally  system  (7)  should  now  be  solved  to 
determine  the  body  shape  l  or  the  tangent  gas  approximation  considered  next  the 
problem  can  be  solved  by  an  exact  method  similar  to  the  one  in  the  incompressible 
case  ( 4  ]. 


5.  Tangevi  Gas  Approximation 


The  tangent  gas  approximation  is  given  by  (see  [21  ]) 


From  (10)  we  obtain 


(24) 


(25) 


where  the  subscript  a  denotes  a  suitable  reference  point  With  the  constant  e  in  (8) 
taken  as 


<  ---  I  'll.,. 


(26) 


we  obtain  from  (S| 


A ( v ) =  1  (27) 

Then  for  subsonic  Don  the  system  (7)  becomes  the  Cauchy  Riemann  equations 


316 


DARIPA  AND  SIROVICH 


Equations  (28)  are  exact  for  both  the  tangent  gas  and  also  for  incompressible  flow 
(A/  =  0).  Henceforth,  we  take  the  tangency  point  (/>„,  !//>„)  to  be  the  free  stream, 

Pu=  P ,  =1'  Pu  =  P  .  =  1  (29) 


''  I—  /f,cothv 


.  (30) 


With  this  selection  the  following  relations  hold  [21  ] 

q  =  sinh  v*  cosechfv*  —  v ),  /f  =  tanh(v*  —  v),  t-„  =  - 

where  the  constant  v*  is  given  by 

m 1 

From  (6)  it  is  seen  that  v,  =0  and  at  a  stagnation  point  (denoted  by  zero  sub¬ 
script) 

s 


It  follows  from  (28)  that 


*'<»=  — x,  c.„  = 


■=  -v  +  iO 


I  +  /*< 


(32) 


(33) 


is  an  analytic  function  of  vr  and  hence  of  n.  A  convenient  representation  of  r(rr)  is 
given  by  (sec  also  [8]), 


exp(r(«r))  =  ( I  —  a  ')  ■'(<'  "'-n  1 )  'exp  (  V  c„< r  "j,  (34) 

where  <)  =  0,/n,  0,  the  trailing  edge  angle.  The  complex  constants  c„  are  denoted  by 

c„  =  A  „  +  (35) 

Note  that  (34)  contains  the  Kutta  condition.  Two  Schwarz  ChristofTel  factors 
appear  in  (33)  because  of  the  discontinuity  in  0  at  the  two  stagnation  points.  On 
the  unit  circle,  (33)  reduces  to 


where 


cxp(T(e",))  =  0'(a)e""’'cxp 


<#(*) 


2  sin  - 


|2(sin  a0  + sin(a  -  a(1))|  1 


'/(*) 


—  <5 )( n  —  y. )  + 


7T (  ( r  —  a,)  +  x,,. 


(36) 


(37) 

(38) 
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U(a-as)  in  (38)  is  the  unit  step  function.  The  tangent  angle  ()„ 
related  to  0  by 

at  the  body  is 

0(a)  =  0B(a)  —  n  —  n(J(a  -  at). 

(39) 

Separation  of  (36)  into  real  and  imaginary  parts  leads  to 

r 

v(a)=  Y.  cos  na  +  B„  sin  na) 

n  =  0 

(40) 

and 

ex. 

8(a)  =  Y  (S„  cos  na  —  A„  sin  na)  +  n  +  aQ 

(41) 

n  ~  0 

where 

v(a)  =  — v(a)  —  In  G(a) 

142) 

and 

t7(a)  =  «fl(a)-^  (1  -<5)(7t-a)-fa  +  ^j. 

(43) 

Notice  from  (34)  the  upstream  flow  direction  0  r  is  related  to  B() 

by 

0 ,  =  Ba  +  n  +  2a0, 

(44) 

where  B{)  is  given  by 

1  C2n 

5„  =  ~  8(a)da-n-aa. 

in  Jo 

(45) 

The  free  stream  condition  (q,,  =  1 )  is  given  by 

/4o  =  0.  (46) 

The  condition  for  closure  of  the  airfoil  is  related  to  the  leading  terms  of  the  series 

(40)  and  (41 )  by  (see  [21  ]) 

A,  =  (1  -<5)-(l -/?«,  )2sin2a0,  (47) 

B,  =  (1  -/?r.)sin  2a0.  (48) 


6.  Behaviour  at  Stagnation  Points 

It  is  both  interesting  and  useful  to  study  the  behavior  of  speed  at  the  stagnation 
points.  From  (30)  and  (42)  we  obtain 


</4 - 2l^—~  for  qx~  0. 

1  +  P G 


T-R^jr’urrRwjnw7XXTrj»n*7\*  vryxTv.rvxrjw  -wuwmr^i 
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From  (37) 


I  fa,'(2cosa0)  for  a~0 

G(ae)  (|a,  -  a|(2  cos  an)*  for  a~a,. 


From  (49)  and  (50) 


A,  a1'  for  a~0 

A\|a, —  a|  for  a~a,, 


e  °’(2cosa0),  A,  =  -  --  e  v~,(* '  J|l(2  cos  a0)\ 


If  a,  and  a,  are  close  to  a^O,  then  from  (51 )  we  obtain 


ln(^((aj)/<7,(a,)) 


In(a3/a, ) 


From  (51 ),  (52),  and  (53)  we  obtain 


f  !«  »  0>_  -In  '-^—aan  A,  I  1 


where  a,  is  close  to  a,  and  a,  is  close  to  zero. 


7.  Mkthod  oi  Soi  tn ion 

The  speed  distribution  </,(.v)  is  usually  given  at  a  finite  number  of  points  .v,.  /  = 
0,  1,  2,...,  in  the  interval  0  ^  v  ^  I .  From  this  the  integral  in  Fq.  (IS)  can  be  evaluated 
to  obtain  (?(.?,)  as  a  function  of  r/J.v,).  Next  the  circulation  /  is  computed  from  the 
relation 

/  =  2(7(.v sss  v, )  Q{s  =  I ). 

Equation  (22)  is  then  solved  to  obtain  a(1  The  second  equation  of  (20)  is  next  used 
to  calculate  the  value  of  the  constant  </.  In  general  a,,  and  a  so  obtained  do  not 
satisfy  the  first  equation  of  (20)  exactly  because  a,,  is  calculated  numerically.  If 


(7(.v  I)  -  (;(v  a,)..  )  </,  <h 


m 


/•  ."*»  ,*•  ,S  ’»  ,%  /•  _*•  ,*»  /•  ^  /•  «  •  ,*•  w  m  M  »  *  •  *  it  **  •  *  "  J»  *  •  "  ■  *Ji  "J* 

V  .<r<-  .’N-  s'f.-A  •-  VaV.sV._-.  .V.V.VvVv.s.y  A  A  A  A 
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differs  slightly  from  Q(s  =  sf)  -4na  sin  a0  (see  Eq.  (21 )),  the  speed  q,(s)  is  modified 
by  a  constant  factor  over  the  interval  s/<s<  1  to  adjust  the  above  integral  lo  this 
value. 

The  values  of  Q(ol,)  at  N  (a  power  of  2)  equally  spaced  points  on  (he  unit  circle, 
a^lnj/N,  y  =  0,  1,...,  N,  are  calculated  using  Eq.  (19).  The  value  of  speed  <y,(a()  at 
the  grid  points  a-  are  now  easily  obtained  by  interpolation  since  </,(.«,)  is  already 
known  as  a  function  of  Qfs,). 

The  approximate  value  of  the  trailing  edge  angle  S  is  obtained  from  Eq.  (53) 
v,(a)  is  then  obtained  from  Eqs.  (42),  (37),  and  (30)  and  its  value  at  a  stagnation 
point  is  calculated  from  (54)  and  (55).  If  v,(a)  satisfies  the  constraints  (46),  (47), 
and  (48)  then  the  conjugate  function  3(a)  is  calculated  from  (41)  using  the  fast 
founer  transform.  In  ^ase  v,(a)  does  not  satisfy  the  constraints,  the  prescribed  speed 
distribution  must  be  modified.  This  is  discussed  in  the  next  section 

The  value  of  the  constant  B0  in  (41)  which  is  also  needed  to  calculate  3(a)  is 
obtained  from  (44)  by  setting  the  free  stream  direction  0,  to  /.cro.  The  tangent 
angle  0H  at  the  body  is  now  obtained  from  3(a)  using  the  relation  (43).  The  body 
coordinates  are  then  calculated  from 


ds 

■  (  cos  0  y. )  d*} f, 
dx 


( 56a ) 


f 5  ds 

r(a)  =  —  sin(>*(a) 

dx 


dx. 


( 56b ) 


where  ds/dx  is  given  by 


ds  .  |sin  a„  +  sin(a  -  a,,)! 

—  =  2  a - 

dx  </,  (a) 


(57) 


The  value  of  (57)  at  a  stagnation  point  (a  =  0,  a  =.  a,)  is  given  by  (see  Eq.  (51 )) 


ds 

dx 


hi 


cos  a,, 


2  a 


cos  a„ 

K: 


for  a  =  0 


for  a  =  a,. 


(58) 


Instead  of  calculating  <)  from  Eq.  (53)  as  was  done  above,  one  can  prescribe  <5 
because  the  constraints  (46),  (47),  and  (48)  depend  on  <3.  Modification  of  the  speed 
distribution  subject  to  these  constraints  will  automatically  satisfy  the  Eq  (53) 
because  this  equation  is  valid  if  the  speed  distribution  is  consistent  with  those  con¬ 
straints.  In  cither  case  if  f,(a)  obtained  from  a  given  speed  distribution  does  not 
satisfy  the  constraints  then  the  prescribed  speed  distribution  must  be  modified 
according  to  the  method  discussed  in  the  next  section 

Even  though  the  above  method  is  exact  theoretically,  there  are  numerical  sources 
of  error.  These  errors  depend  on  the  kind  of  interpolation  and  integration  scheme. 


■  »  ..t  l^i  M  I 


U|  l.l  ».  I 
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the  number  of  data  points,  the  number  of  grid  points,  the  evaluation  of  a„  from  (22) 
and  the  use  of  the  approximate  expressions  (53),  (54),  and  (55)  to  calculate  trailing 
edge  angle  <5,  v,(a  =0),  and  v,(a  =  a,),  respectively. 

In  view  of  the  simplicity  of  the  procedure  no  attempt  was  made  to  incorporate 
highly  accurate  computations.  Simpson's  rule  with  evenly  spaced  grid  points  and 
trapezoidal  rule  with  unevenly  spaced  data  points  were  used  for  integration.  The 
interpolation  scheme  used  was  linear.  The  speed  </,(.v)  was  prescribed  at  129 
unevenly  spaced  data  points  on  0$.!^  I  and  the  number  of  grid  points  on  the  unit 
circle  was  taken  to  be  128.  a„  was  obtained  within  an  accuracy  of  10  6  by  solving 
Eq.  (22)  by  regular  falsi  method  and  the  trailing  edge  angle  d  used  was  calculated 
by  using  the  approximate  relation  (55). 

The  program  was  run  on  an  IBM  3081  in  single  precision  and  the  computation 
time  was  about  a  half  second  in  most  cases. 


8.  Modii  ic'A  i  ion  oi  Sim  i  d  Distribution 

Constraints  (46),  (47),  and  (48)  must  be  satisfied  by  the  prescribed  speed  dis¬ 
tribution  to  find  a  closed  body  solution.  Therefore  in  general  any  arbitrary  spiced 
distribution  must  be  modified  subject  to  these  constraints.  These  constraints  can  be 
written  in  terms  of  surface  values  v',(x)  given  bv 

—  j  v,(  7 )  g. ( 7 )  (h  =  j  1 .  2.  3.  ( 59 ) 

r, 

where  ,tf;(a)  and  l\  are  given  by  (see  (46),  (47),  and  (48)) 
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where  yk,  A  =  I,  2,  3,  are  constants  to  be  determined  and  /*(a),  A  =  1,  2,  3.  are 
suitable  correction  terms.  The  correction  terms  can  be  set  to  zero  outside  a  specified 
interval  (a,,  a,)  leaving  speed  distribution  same  as  the  prescribed  one  outside  this 
interval.  This  is  extremely  useful  when  designing  an  airfoil  where  in  general  no 
modification  of  speed  distribution  over  the  suction  side  is  desired.  The  speed  dis¬ 
tribution  can  be  modified  in  various  ways  depending  on  the  choice  of  functions 
/*(ot),  A  =  I,  2,  3.  and  the  correction  interval  ( 3t , ,  a,). 

Substituting  (62)  in  (59)  one  obtains 


Z  7*",‘  =hr  /=  1.  2.  3.  (63) 

A  1 

where 

'6*  =  I  /i<2)  /*( 3t )  g,(x)  ih.  (64) 

Ml 

hl-nPl-\  vjot)  ,i,',(a)Jx  =  nl\  - j  vji)  g,(a)  </x.  (65) 

•  'i 

Constants  yk ,  A  =  1,  2,  3,  are  obtained  by  inverting  (63)  and  the  corrected  0,(3!)  is 
then  obtained  from  (62).  The  corrected  speed  distribution  is  then  obtained  from 
(42)  and  (30)  and  the  body  is  found  from  (55). 

The  matrix  in  (65)  must  be  positive  definite  to  be  able  to  invert  (h.3)  which 
restricts  the  choice  of /*(a).  a,,  and  ot:.  These  should  be  carefully  selected  so  that 
the  correction  to  a  prescribed  speed  distribution  is  minimum.  This  can  be  done  in 
the  same  spirit  as  in  Strand  [22]  and  Arlinger  [4], 

For  our  purpose  we  choose  a^O.  a:  =  2 n  and  J k(r)  -  .gjx).  A  =  1,2.  3.  This 
modifies  the  speed  distribution  over  the  whole  interval.  In  this  case  (63)  gives 


7 1 


and  hence  (62)  becomes 


I 

v,(a)  -  v ,„( T )  -)-  —  (fr|  +  2h:  cos  a  +  2 />,  sin  a). 


(66) 


(67) 


9. 


Ristii  is 


A  basic  test  of  the  inverse  method  is  the  recovery  of  a  known  airfoil  from  its 
pressure  distribution.  Figures  2  and  3  provide  a  verification  of  the  method  within 
the  tangent  gas  approximation.  Here  a  pressure  distribution  is  computed  in  tangent 
gas  approximation  over  a  NACA  4412  airfoil  (see  [21  ] )  Speed  distribution  is  com- 


i 

i 


« 

I 


« 

y 

\ 


322 


DARIPA  AND  SIROVICH 


»  ♦  ♦  »  PRESSURE  DISTRIBUTION  ON  DESIGNED  AIRJ  OIL 


Fig.  2.  Comparison  of  pressure  distributions  on  NACA  4412  airfoil  (input)  and  on  designed  airfoil 
from  tangent  gas  solution  at  free  stream  Mach  number  0.7  and  zero  angle  of  attack 


puted  from  this  pressure  distribution  using  Kqs.  (30).  Then  the  airfoil  is  designed 
from  this  speed  distribution  by  the  method  discussed  in  Section  7.  In  Fig.  2  we  show 
the  pressure  distribution  over  a  NACA  4412  airfoil  as  calculated  by  the  tangent  gas 
and  compare  it  with  the  pressure  over  the  designed  airfoil.  The  error  is  less  than 
<2(10  ').  Figure  3  compares  the  given  airfoil  with  the  designed  airfoil.  The  error  is 
less  than  0(10  4).  (The  origin  of  these  errors  is  numerical  and  was  discussed  in 
Section  7.) 

Figure  4  shows  a  pressure  distribution  which  did  not  satisfy  the  constraints  (46), 
(47),  and  (48).  The  pressure  distribution  which  results  from  the  correction 
according  to  (67)  is  shown  in  the  same  figure.  The  resulting  body  along  with  its 
design  and  analysis  pressure  is  shown  in  Fig.  5. 


A 


Ik;  3.  Comparison  of  NACA44I2  airfoil  ami  the  airfoil  designed  hv  tangent  gas  approximation 
from  input  pressure  distribution  of  l  ig  2 
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Fig.  4.  Target  and  concerted  pressure  distribution. 
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Fig  5.  Designed  airfoil  from  input  pressure  distribution  of  Fig  4  and  pressure  distribution  on 
designed  airfoil  from  design  and  direct  analysis  b\  tangent  gas  method 
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Fl(i.  7.  Comparison  of  NACA00I2  and  the  airfoil  designed  by  tangent  gas  approximation  from 
speed  distribution  of  Fig.  6. 
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This  speed  distribution  is  used  to  design  the  airfoil  by  the  method  mentioned 
before.  Figure  6  shows  the  speed  distribution  on  a  NACA0012  airfoil  as  calculated 
from  the  Euler  equations  at  free  stream  Mach  number  0.6  and  zero  angle  of  attack. 
Figure  7  compares  the  NACA0012  airfoil  and  the  designed  airfoil  in  the  tangent  gas 
approximation.  It  is  seen  that  the  airfoil  is  almost  exactly  recovered  along  with  the 
zero  angle  of  attack.  The  pointwise  error  is  less  than  3  %  and  this  only  occurs  in  a 
small  neighborhood  of  the  leading  edge.  Figure  8  compares  the  Euler  pressure  with 
the  pressure  over  the  designed  airfoil.  Again  the  comparison  is  excellent  except  near 
the  leading  edge  where  the  error  in  Cp  is  0(10  2).  It  is  to  be  emphasized  that  this 
error  occurs  as  a  result  of  using  the  tangent  gas  approximation  and  is  in  no  way 
numerical.  We  believe  on  the  basis  of  this  discussion  that  this  recommends  the  use 
of  the  method  presented  here  for  airfoil  design  especially  since  it  is  computationally 
very  efficient. 

In  the  next  example  we  push  the  method  beyond  its  limits  by  considering  a 
supercritical  case.  We  show  in  Fig.  9  the  Euler  speed  distribution  over  a  NACA0012 
airfoil  at  free  stream  Mach  number  0.5  and  angle  of  attack  5  .  F'igure  10  shows  that 
we  recover  the  correct  angle  of  attack  and  the  airfoil  except  over  a  small  region 
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H(;.  8.  Comparison  of  Fuler  pressure  distributions  over  NACA0012  airfoil  and  the  airfoil  designed 
from  speed  distribution  of  Fig  <>  at  M  ,  O.ft  and  r  0.0. 
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I  k,  9 

j  ^  5 

near  the  nose  where  the  error  is  within  2%.  and  where  the  flow  is  actually  super¬ 
critical.  Again  the  error  is  nonnumcricai  and  gives  a  measure  of  the  deviation  of  the 
tangent  gas  approximation  from  the  exact  Euler  result.  In  Pig.  1 1  we  compare  the 
Euler  solution  over  the  NACA00I2  airfoil  and  designed  airfoil.  Note  that  the 
agreement  near  the  leading  edge  is  not  as  good  as  elsewhere  because  the  designed 
airfoil  suffers  maximum  deviation  from  the  NACA0012  airfoil  near  the  leading  edge 
Finally  a  useful  application  of  our  approximate  method  is  to  provide  a  starting 
airfoil  in  a  design  procedure  in  which  the  Euler  equations  are  used  directly  to  give 
the  corrected  pressure  distribution.  At  successive  stages  the  pressure  distribution  is 
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I  k,  1 1  Comparison  of  l-ulcr  pressure  distributions  over  NACA00I2  airfoil  and  the  airfoil  designed 
from  speed  distribution  of  l  ie  9  at  M  ,  -0.5  and  i  5 


modified  to  meet  the  design  criteria  and  the  inverse  method  reapplied  and  so  forth. 
The  interactive  iteration  should  go  quite  quickly  for  subcritical  flows.  However,  the 
raiuc  of  such  a  procedure  remains  uncertain  for  supercritical  flows. 
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FROM  GIVEN  PRESSURE  DISTRIBUTION  IN 
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Abstract 


A  new  inverse  method  for  aerodynamic  design  of 
airfoils  is  presented  for  subcritical  flows.  The  pressure 
distribution  in  this  method  can  be  prescribed  in  a  natural 
way.  i.E.,  as  a  function  of  arclensth  of  the  as  yet 
unknown  body.  This  inverse  problem  is  shown  to  be 
mathematically  equivalent  to  solving  only  one  nonlinear 
boundary  value  problem  subject  to  known  Dirich'et  data 
on  the  boundary  The  solution  to  this  problem  deter¬ 
mines  the  airfoii.  free  stream  Mach  number  M«  and  the 
upstream  liow  direction  0,  The  existence  of  a  solution 
to  •  give-  pre-sire  distribution  is  discussed.  The 
method  :s  ru-v  :r:  mplemer.t  and  extremely  efficient. 
W*  vese-:  :  -ene-  >r  results  for  which  comparisons  are 
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boundary  layer  separation,  wave  drag  and  lift  Success 
with  this  method  in  designing  cost  efficient  airfoils  has 
made  this  approach  useful  in  airfoil  design. 

The  mathematical  theory  behind  solving  the  inverse 
problem  is  somewhat  more  difficult  than  the  correspond¬ 
ing  analysis  problem  for  the  following  reasons:  a)  diffi¬ 
culty  with  formulation  of  proper  differential  equations: 
b)  lack  of  existence  of  a  solution  for  arbitrarily 
prescribed  pressure  data;  c)  the  difficulty  in  imposing 
closure  constraints  (a  specified  gap  at  the  trailing  edge  of 
the  designed  airfoil).  For  incompressible  flows  ihe 
theory  of  harmonic  functions1-2  makes  the  above  issues 
easily  tractable  due  to  the  linearity  of  the  problem  Flerc 
the  existence  and  the  closure  constraints  can  easily  be 
established  The  specified  pressure  distribution  can  be 
modified  apriori  to  satisfy  these  constraints  Most  of  the 
solution  techniques  tor  this  case  have  been  based  on  ana¬ 
lytic  function  theory;'9  The  tangent  gas  approximation 
makes  the  inverse  problem  for  subsonic  flows  similar  to 
the  incompressible  case10,11 

Nonlinearity  in  compressible  flows  makes  this  prob¬ 
lem  difficult.  Most  of  the  inverse  methods  rely  on  either 
a  Dmchlct  or  Neumann  formulation  depending  on  the 
choice  of  the  dependent  variable  12'!S  and  usually 
involves  solving  a  series  of  nonlinear  elliptic  problems 
An  excellent  account  of  these  methods  can  be  found  in 
Sloof1,9  Formulation  of  the  closure  constraints  and  the 
existence  of  a  solution  for  such  flows,  similar  to  the  case 
of  incompressible  Hows,  may  prove  useful  in  developing 
new  numerical  techniques  The  constraint  nece-siated  by 
the  existence  of  a  solution  in  the  subsonic  case  was  .-sia 
hlishcd  :n  Darip.i~.u  It  is  shown  '.here  that  because  >0  -he 
nonlinearity  ol  the  constraint,  the  existence  o:  a  solution 
.an  not  be  established  ipri  -ri  from  the  pre-cr  -cd  .  -.  s 
sure  distribution  However  u  suggests  ;ha:  the  a-  -  ■  : 
formulation  o!  ihe  problem  nii-l  base  at  lea-.'  me  ':■■■■ 
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parameter  to  make  the  problem  solvable  in  general. 
Thts  tree  parameter  will  be  determined  by  the  solution. 
Such  an  attempt  has  been  made  in  two  dimensions'-  for 
transonic  flows. 

The  case  of  transonic  flow  is  more  difficult  than  the 
other  cases  owing  to  the  mixed  elliptic-hyperbolic  nature 
of  the  transonic  flow-  equations.  The  mathematically 
elegant  method  of  complex  characteristics,  was  success¬ 
fully  used  by  Garabedian-1'"4  to  generate  supercritical 
a tr foils  In  this  method  the  boundary  is  unknown  and 
iteration  on  the  boundary  is  used  to  generate  the  airfoil. 
There  are  also  cost  effective  methods  based  on  the  ficti¬ 
tious  gas  concept  to  generate  supercritical  airfoils,-5  An 
excellent  review  on  the  design  of  supercritical  airfoils 
and  wings  can  be  found  in  Sobieczky-.6 

v\  e  discuss  our  method  in  this  paper  for  subsonic 
[lows  only  In  our  formulation,  the  equations  of  motion 
are  cast  into  one  boundary  value  problem  in  the  poten¬ 
tial  olar.e  To  avoid  dealing  with  the  infinite  potential 
plane,  we  map  our  potential  plane  into  the  interior  of  a 
unit  circle  and  solve  our  equations  there.  To  render  the 
problem  solvable  we  choose  the  Mach  number  distribu¬ 
tion.  computed  from  the  prescribed  pressure  distribu¬ 
tion,  as  the  boundary  data.  By  doing  this  we  have  at  our 
disposal  the  free  stream  Mach  number  as  free  parameter 
which  ;s  determined  as  part  of  the  solution  of  the  boun¬ 
dary  value  problem.  The  solution  also  determines  the 
shape  of  a  profile.  This  free  stream  Mach  number  and 
the  input  Mach  number  distribution  together  determine 
the  pressure  distribution  over  the  designed  profile.  In 
generai  the  prescribed  and  computed  pressure  distribu¬ 
tions  wtii  not  be  same  except  in  cases  where  there  exists 
a  solution  :o  the  prescribed  pressure  distributions. 

Our  method  requires  solving  only  one  nonlinear 
boundary  value  problem  as  opposed  to  solving  a 
sequence  of  such  problems  In  addition  a  profile  is 
always  generated  by  our  method  The  designed  airfoil 
however  may  have  a  gap  at  the  trailing  edge.  To  be  able 
to  design  an  airfoil  with  any  prescribed  gap  we  need  to 
do  further  work  that  :s  currently  in  progress  An  obvi¬ 
ous  approach  is  to  use  this  method  in  some  iterative 
mode 

II.  Prohlem  Formulation 

The  eueil:o-  ol  motion  arc 

V  •  pq  i  <i  V  «•  q"  0  .  p  -  py  tla.b.e) 


The  variables  are  normalized  by  their  some  values  and 
linear  dimensions  by  some  appropriate  linear  dimension 
Here  p  is  the  pressure,  p  is  the  density  and  q  is  the 
speed.  These  equations  imply  the  existence  of  a  stream 
function  vli  and  a  potential  6  given  by 

pq  =  Vx(ilzF)  :  q  =  V4>.  (2a. b) 

where  £  denotes  a  unit  vector  perpendicular  to  the  plane 
of  the  motion.  The  above  equations  can  alternatively  be 
put  in  the  form 

9*  -  (K(Mirl»t  =  0  :  8*  +  (K(M)iv*  =  0 

f  3a . b  J 

Here  K  and  v  are  functions  of  Mach  number  M  only-' 
and  8  is  the  flow  direction,  v  is  known  as  the  Prandtl- 
Mever  function.  The  body  maps  onto  a  slit  in  the  poten¬ 
tial  (w  =  <f>  +  i  xlz)  plane  as  shown  in  Fig.  1. 

Differentiation  and  elimination  reduces  the  system 
(3a. b)  of  first  order  partial  differentia]  equations  (PDEs) 
to  an  equivalent  second  order  PDE  in  v  only. 

(K-'tMlv^  +  (K(M)w4)4  -  0. 

A  little  algebra  shows  that  this  equation  can  alternatively 
be  written  in  the  form-' 

=  (1  ~  K:)vi6  +  flMfvj,'-  -  K:f(M/vj,:  ( 53 
where  f(M)  is  a  function  of  Mach  number:' 

The  appropriate  boundary  condition  is  seen  to  be 
(see  Fig.  I) 

for  viz  =  0 
for  di  =  0 

where  (tb.ili)  =  (d)A  =  0,0)  corresponds  to  the  front 
stagnation  point  and  (6.Jz)  =  i  dig, 01  and 

(d>,Jz)  =  (ctjy.O)  correspond  to  the  upper  and  lower  side 
of  the  rear  stagnation  point  respectively  (see  Fig  1) 
Bernoulli's  law  gives  q  =  q(Cp)  and  M  =  M(Cpi  -s  -u 
that  are  used  to  determine  q(s)  and  M(s)  from  the  input 
pressure  distribution  Cp(s);  0  5  s  <  1.  Here  s 
parametrizes  the  arclength.  The  equation  (2b)  implies 

<b  =  J!q(s)dsi  <"' 

which  together  win  known  q(s)  and  M(s),0  £  >  s  I 
determines  the  boundary  condition  (6)  and  the  boundary 
(the  slit)  in  the  potential  plane 


M=.M(6) 
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IV.  Results 


The  solution  of  the  boundary  value  problem  (5)  and 
(6)  determines  v0  on  the  slit  that  is  subsequently  used  in 
(3a)  to  compute  94  on  the  slit.  Integrating  this  and 
using  (7)  determine  the  flow  angle  as  a  function  of 
arclength  on  the  body  and  hence  the  profile  is  known. 

To  avoid  dealing  with  infinity  in  the  (4>  -  il/i  plane 
we  map  this  plane  into  the  interior  of  a  unit  circle  such 
that  the  body  maps  onto  the  unit  circle.  We  carry  out  our 
calculations  in  this  circle  plane  }' 

The  solution  of  this  boundary  value  problem  also 
determines  the  upstream  Mach  number.  Pressure  on  the 
designed  body  can  be  computed  using  the  computed  Nl, 
and  the  input  Mach  number  distribution.  The  flow 
direction  at  infinity  with  respect  to  the  body  can  be  cal¬ 
culated  by  integrating  either  f  3 a i  or  (3b)  from  the  siit  to 
infinity  in  the  potential  plane.  This  integration  process 
becomes  much  easier  in  the  mapped  circle  plane.  See 
Daripa  for  details 

III.  Method  of  Solution 

As  mentioned  earlier,  the  boundary  and  the  boun¬ 
dary  data  in  the  potential  plane  and  hence  in  the  mapped 
circle  plane  are  known  from  the  input  pressure  distribu¬ 
tion.  The  equation  (5)  is  solved  numerically  inside  the 
unit  circle  subject  to  the  known  Di  rich  let  data.  In  solv¬ 
ing  equation  if)  numerically,  a  general  linear  Helmholtz 
equation  solver  j0  is  used  iteratively.  An  initial  guess  of 
the  flow-field  inside  the  unit  circle  determines  the  right 
hand  side  of  equation  (5)  within  the  unit  circle  The 
solver  then  updates  the  values  of  v  inside  the  unit  circle 
by  solving  the  linear  Poisson  equation  and  the  process  is 
repeated  until  a  given  convergence  criterion  i-.  met 

The  solution  is  considered  to  have  converged  if  the 
difference  in  the  maximum  value  of  Prandtl-Mever  iunc- 

tion  v  inside  the  unit  circle  between  two  successive 
iterations  is  less  than  5X10-6  .  We  never  needed  .  re 
than  six  iterations  in  any  ol  our  calculations  when  the 
initial  guess  was  taken  to  be  that  of  uniform  flowtield 
The  solution  values  of  v  inside  the  unit  circle  are  then 
used  to  compute  the  estimate  of  the  normal  derivative  on 
the  unit  circle  This  estimate  is  then  used  to  compute  the 
body  angle  through  the  use  of  equation  (3a)  and  the 
mapping  function  Similarly  the  (low  direction  at  infinity 
is  calculated  by  integrating  equation  (3ai  along  an 
appropriate  ray  in  the  circle  plane.  The  solution  ot  the 
nonlinear  elliptic  equation  (f)  also  determines  the 
upstream  Mach  number 


To  validate  our  method  we  present  a  series  at 
results  which  recovers  a  known  airtoil  from  its  pressure- 
distribution.  We  generate  pressure  distributions  over  a 
scries  of  closed  airfoils  at  a  given  tree  stream  Mach 
number  M*  and  angle  of  attack  9,  by  using  an  Euler 
code  (  flo52s  written  by  A.  Jameson.  E  Turkel  and  M 
Salas  ).  This  pressure  distribution  is  then  used  in  our 
method  to  generate  the  airfoil.  Numerical  sources  of 
error  in  practice  may  introduce  some  error  and  so  we 
monitor  the  following  as  a  measure  ot  accuracy  ot  our 
method:  E(M)=  |  .  E(8,)  =  |0,  -  6i !  . 

E(9b)  =  max  |93(s,)  —  9§(s,)  |  .  veap  =  |y(s  =  0)  -  yi  s  =  1 )  i 
and  xgap  =  !.x(s  =  0) -x(s=  1)  |  .  where  xgap  and  yeap 
(normalized  by  chord  length)  measure  the  gap  at  the 
trailing  edge  of  the  airfoil  and  90  refers  to  the  body 
angle.  In  the  above  a  superscript  refers  to  the  computed 
value.  Here  we  present  a  few-  results  for  cases  where  the 
trailing  edge  is  not  a  stagnation  point.  The  case  with  rear 
stagnation  point  has  been  excluded  here  because  some 
difficulties  was  encountered  in  removing  the  singularity 
at  the  trailing  edge.-7  This  case  will  be  taken  up  m 
future. 

Fig.  2  shows  the  Euler  pressure  distribution  over  a 
12%  thick  Kutta  airfoil  at  M*  =  0.5  and  angle  of  attack 
=  2.0  degrees.  The  application  of  the  above  method 
then  generates  the  body  and  also  gives  the  computed 
values  of  the  free  stream  Mach  number.  M i .  anti  the 
angle  of  attack.  9£  The  number  of  iterations  ( -ee  t;  III) 
required  to  converge  to  the  solution  using  the  linear 
elliptic  solver  were  only  five  in  this  case.  In  Fig  a  vv  e 
compare  the  designed  airfoil  with  the  exact  airtoii  and 
find  that  the  agreement  is  excellent.  We  rind  the  com¬ 
puted  free  stream  Mach  number  MS  -  0  50005  and  -he 
angle  of  attack  9£  —  2  05  degrees  The  value-*  ot  the 
error  diagnostics  in  this  case  are  E(Ml  =  0  OOOOn  rA, 
)  =  0.00087.  E(90)  =  0.009,  xgap  =  0.000-13  and  ygap 
=  0.00010.  The  cap  at  the  trailing  edge  is  within  0  01% 
of  the  thickness  of  the  airfoil  which  is  negligibly  '.mall 
A  more  important  quantity  is  the  body  angle  which  is 
likelv  to  suffer  maximum  error  near  the  leading  edge 
since  the  body  angle  is  a  rapidly  varying  (unction  of 
arclength  there.  Fig.  -!  compares  the  computed  and  exact 
values  of  body  ancle  as  a  function  of  arclength  I  i  -  s 
compares  the  same  in  the  leading  edge  region  Here 
again  we  hnd  the  error  in  the  body  angle  w  ()  I" 
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every  w  her ■•  except  near  the  leading  edge  where  it  is 
maximum  and  is  given  by  E  iSp)  =  0  009  In  general 
closure  can  not  be  expected  unless  explicitly  imposed 
either  by  some  iterative  or  exact  method  i\Ve  will 
report  on  the  design  of  closed  airfoil  in  the  near  future.) 
The  source  o;  these  errors  is  numerical  as  has  been  dix- 


Nex:  iv;  mow  a  more  vigorous  test  case  for  gen¬ 
erating  an  airfoil  from  an  arbitrary  pressure  distribution 
Seen  an  arbitrary  pressure  distribution  (solid  iine)  is 
shown  in  F'c  6  The  asymmetric  pressure  distribution 
in  Fig  6  is  a  minor  variation  of  the  pressure  distribution 

;p  rig.  I  and  as  we  shall  'C-e  there  no  airfoil 
corresponding  to  this  input  pressure  distribution  The 
resulting  airfoil  designed  with  our  inverse  method  is 
shown  in  rig  7.  This  airfoil  has  a  t  in  tie  gap  at  the  trail¬ 
ing  edge.  We  find  the  computed  Mach  number  M£  = 
0  sjtu  and  angie  of  attack  Si  =  2  degrees  The  pressure 
distribution  on  the  designed  airfoil  obtained  by  our 
inverse  cooe  a:  this  Mach  number  and  angle  of  attack  is 
also  show  n  in  Fig.  6  (  a-  sign)  which  is  different  from  the 
input  pressure  distribution  (solid  line).  This  is  due  to  the 
tact  that  there  is  no  airfoil  associated  with  this  input 
pressure  distribution .  In  order  to  have  a  specified  gap  at 
the  trailing  edge,  this  method  can.  be  used  in  an  iterative 


V.  Discussions  and  Conclusions 

Vs  e  ha'.e  developed  a  fast  approach  to  solve  inverse 
problem  for  -abcritica!  flows  We  have  shown  that  solu¬ 
tion  of  the  -verse  problem  requires  solving  only  one 
nonlinear  boundary  value  problem  in  an  appropriate 
plane  Mach  number  distribution  is  used  as  the  boun¬ 
dary  vaiues  which  arc  in  turn  used  to  generate  the  atr- 
tiui  i  bis  approach  leaves  the  upstream  Mach  number  as 
j  tree  parameter  which  is  also  determined  bv  the  >olu- 


The  numerical  sources  of  error  may  introduce  a  gap 
ut  the  trailmc  edge  even  in  cases  where  ideally  there 
-hue Id  be  no  cap  Therefore  it  is  necessary  to  mcor- 
•sorate  >ome  -.ird  ot  closure  procedure  m  our  present 
tiv—.  nation  it  is  my  .intention  'hat  this  method  can 
easii'-  be  exiemted  to  deal  with  the  design  ot  a  supercriti¬ 
cal  airtnii  n  in  etticien:  manner  However  such  an 
attempt  has  ro  t  vet  been  nude 
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Fig.  1  Airfoil  in  the  physical  tt-pluno 
and  potential  w-plane. 


’sr.v.v.-i  -■ 


.4  .6  .3 

Ardength 


j'pper  sid»- 


Fig.  2  Pressure  distribution  over  a  12%  thick  Kutta  airfoil 
from  Euler  solution  (fio52s)  at  M«  =  0.5  and  a  =  2°  degrees. 
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.  Body  angle  of  the  Kutta  airfoil 

+  +■  Body  angle  of  the  airfoil  designed  from  pressure 
distribution  of  fig.  2. 


Fig.  4  Comparison  of  the  body  angle  of  the  designed  airfoil 
and  of  the  Kutta  airfoil  as  a  function  of  ardcngth  of  the  air¬ 
foil.  The  input  pressure  distribution  of  the  designed  airfoil  is 
shown  in  Fig.  2  . 
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Fig.  3  Comparison  of  the  12%  thick  Kutta  airfoil  and  the  air¬ 
foil  designed  from  the  Euler  pressure  distribution  of  Fig.  2  . 
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.  Body  angle  of  the  Kutta  airfoil  in  the  leading  edge 
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Fig  5  Comparison  o(  the  bods  angl-  of  the  designed  airfoil 
and  of  the  Kutta  airfoil  as  a  function  of  ardrngih  of  the  air¬ 
foil  in  the  leading  edge  region  The  input  pressure  distribu¬ 
tion  of  the  designed  airfoil  is  shown  m  1 1 h  2 
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